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ABSTRACT 

An  extensive  study  of  an  existing  digital  computer  program 
for  determining  parameters  required  to  represent  SES  100B  trials  data 
by  a  generalized  gamma  function  is  carried  out.  The  purpose  of  the 
study  is  to  resolve  the  difficulties  encountered  in  earlier  studies. 

In  these  earlier  studies,  approximately  5  percent  of  the  cases  analyzed 
produced  unrealistic  values  for  the  distribution  parameters.  In  order 
to  resolve  this  problem,  data  sensitivity  studies  are  carried  out. 

Based  on  results  of  the  study,  criteria  are  developed  for  determining 
cases  that  should  not  be  processed.  Three  other  techniques  besides 
the  Stacy-Mihram  method  are  investigated.  The  most  successful  alterna¬ 
tive  method  is  the  maximum  likelihood  method,  which  is  incorporated  into 
the  new  digital  computer  program.  A  program  listing  and  an  example 
of  the  generated  output  are  included  in  this  report. 

ADMINISTRATIVE  INFORMATION 

This  work  was  conducted  at  the  David  W.  Taylor  Naval  Ship 
Research  and  Development  Center  and  authorized  by  Naval  Sea  Systems 
Command,  Surface  Effect  Ships  Project  Office,  Program  Element  Number 
63534N.  It  is  identified  as  Work  Unit  Number  1-1506-012. 

INTRODUCTION 

A  computerized  procedure  for  evaluating  the  parameters  of  the 
generalized  gamma  distribution  from  a  set  of  random  data  following  the 
method  proposed  by  Stacy  and  Mihram  (see  Reference  1)  has  been  previously 


developed  and  extensively  exercised  (see  Reference  2  and  3).  During 
these  earlier  studies  it  was  discovered  that  what  appeared  to  be  unreal¬ 
istic  values  for  the  distribution  parameters  were  obtained  in  roughly  5 
percent  of  the  cases  analyzed.  The  initial  purpose  of  the  study,  whose 
results  are  summarized  in  this  report,  was  to  explore  in  detail  a  repre¬ 
sentative  sample  of  cases  resulting  in  unrealistic  values  of  the  param¬ 
eters,  and  if  possible,  to  develop  an  alternate  method  for  processing 
these  cases. 

The  results  of  the  present  investigation  can  be  summarized  as 

follows: 

a.  In  the  many  cases  studied,  it  appears  that  the  predicted 
extreme  design  values  are  most  sensitive  to  experimental 
measurements  and  not  the  numerical  procedure  for  deter¬ 
mining  the  parameters.  This  is  particularly  true  when 

one  of  the  parameters  is  negative. 

% 

b.  All  four  methods  considered  in  this  investigation 
yield  probability  distribution  functions  which  represent 
very  well  the  experimentally  obtained  histogram  even 
though  some  of  the  parameter  values  were  judged  unreal isi tic 
in  previous  studies.  However,  the  Stacy-Mihram  and  maximum 
likelihood  methods  are  by  far  the  most  efficient  and  are 
the  preferred  procedures  based  on  probabilistic  arguments. 

c.  In  some  cases  the  predicted  extreme  design  values  are 
unrealistically  large  whereas  the  significant  values  are 
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realistic.  Based  on  an  analysis  of  over  100  cases,  it 
can  be  concluded  that  cases  having  a  nondimensional  third 
logarithmic  moment  greater  than  -0.5  produce  unrealistic 
extreme  design  values.  Results  of  a  data  sensitivity  study 
have  revealed  that  in  all  these  cases  the  data  may  not  be 
considered  as  a  sample  taken  from  a  steady-state  random 
phenomenon,  and  hence  it  is  recommended  not  to  carry  out 
statistical  analysis  for  these  cases. 
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MATHEMATICAL  FORMULATION 


OVERVIEW 

The  generalized  gamma  density  function  under  consideration  in 
this  report  has  the  form 

crxM'l»1'a,l‘)  (,) 

rc«o 

where  X,  m,  and  c  are  estimation  parameters  and  r(m)  is  the  gamma  function. 
Many  other  distributions  such  as  the  exponential,  gamma,  chi-squared, 
Weibul,  hydrograph,  chi,  truncated  normal,  and  Rayleigh  are  special  cases 
of  the  generalized  garrana  function.  Due  to  the  general  nature  of  f(x;  X, 
m,  c),  it  is  the  preferred  distribution  particularly  when  the  exact  form 
of  the  distribution  function  is  unknown. 

There  are  several  techniques  that  can  be  used  to  estimate  the 
three  parameters.  Four  such  techniques  are  considered  in  this  effort. 
These  are  the  Stacy-Mihram,  the  maximum  likelihood,  the  nonlinear 
least  squares,  and  grid  search  methods.  The  last  two  methods  deter¬ 
mine  the  parameters  by  searching  for  that  set  of  parameters  resulting 
in  a  minimum  error  between  the  experimental  and  theoretical  histograms. 
Since  statistical  inference  cannot  be  applied  in  a  rigorous  manner,  the 
last  two  methods  are  not  recommended.  The  Stacy-Mihram  method 
determines  the  parameters  by  requiring  that  the  first  three  theoretical 
logarithmic  moments  agree  with  the  logarithmic  moments  determined  from 
the  experimental  data.  The  maximum  likelihood  method  determines  param¬ 
eters  such  that  the  joint  probability  density  function  is  a  maximum. 
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A  detailed  derivation  of  the  Stacy-Mihram  method  is  presented 


in  Reference  2  pages  6-8.  A  derivation  of  the  nonlinear  least  squares 
procedure  is  presented  in  Reference  4.  The  grid  search  procedure  merely 
involves  a  systematic  search  of  the  X,  m,  c  space  for  parameters  that 
minimize  the  error  between  the  predicted  distribution  and  the  measured 
distribution.  A  derivation  of  the  equations  required  in  the  maximum 
likelihood  parametric  estimation  method  is  presented  in  the  next  section. 
MAXIMUM  LIKELIHOOD  METHOD 

Briefly,  the  maximum  likelihood  estimation  technique  involves 
finding  the  parameter  (or  parameters)  which  maximizes  the  joint  proba¬ 
bility  density  function  of  a  random  sample  .  XN  from  a 

distribution  with  f(x,  0)  as  its  probability  density  function  and  9  as 
an  element  of  the  parameter  space.  The  joint  probability  density  func¬ 
tion,  referred  to  as  the  likelihood  function,  is  given  by 

u-  ir  k  v  s)  (2) 

Assume  one  has  N.  samples  at  random  variable  x^.  Then  the 
i  i 

maximum  likelihood  function  will  have  the  form 


L‘FKx,-,e) 
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(3) 
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The  natural  logarithm  of  the  likelihood  function  is  used  in 


maximizing  since  it  is  a  maximum  when  the  function  itself  is  a  maximum, 
and  it  proves  more  convenient  to  use.  The  logarithm  is  given  by. 


\«N  L  =  1  M-  (4) 

The  maximum  of  InL  is  obtained  by  searching  for  parameters  9  such  that, 

.  =  0.0  (5) 

where  Np  *  the  number  of  parameters. 

For  the  generalized  gamma  function  given  by  Equation  (1), 
the  likelihood  function  becomes 


ainL 


iN; 

^  J 


(6) 


and  InL  becomes 

In  L  »  N  \  n  c  +  H  \a  ^  -  N  U  T  Cm) 

A.  Ccm-0  %  U*-  -  3  %  h  N(- 

l  1  l  1 


(7) 


In  the  above  equation,  c  is  assumed  positive.  If  the  above  equation  is 
differentiated  with  respect  to  cm,  X  and  m  and  set  to  zero,  the  following 
equations  can  be  derived  (see  Reference  5). 


2  s  L  Cm  M)  /&.  ^  m.  ]  ^ 

C 


Once  a  value  of  c  is  determined  by  solving  Equation  (12), 
values  for  X  and  m  can  be  obtained  via  Equations  (8)  and  (11).  The 
numerical  method  used  to  solve  Equation  (12)  for  c  is  due  to  Wegstein 
(see  Reference  5). 

ANALYSIS 

During  previous  studies  in  which  the  Stacy-Mihram  method  was 
used  to  process  SES  1008  data,  large  design  extreme  values  were 
predicted  in  approximately  5  percent  of  the  cases  although  the  most 
significant  values  agreed  with  the  experimental  measurements.  This 
difficulty  seemed  to  be  associated  with  large  values  of  the  parameter 
m  and  with  negative  values  of  the  parameter  c.  In  order  to  resolve 
this  difficulty,  several  cases  were  examined  in  detail.  The  results 
of  this  investigation  follows. 

The  first  case  considered  in  this  investigation  is  a  set 
of  pitching  motion  data  measured  during  the  SES  100B  full  scale  trials 
program  at  30  kts  in  a  seastate  of  3.  This  particular  case,  identified 
as  ART  3E  30,  is  presented  in  Reference  3.  A  comparison  of  the 
estimated  probability  density  function  with  the  experimental  histogram  is 
presented  in  Figure  1.  As  is  demonstrated  by  the  results  in  the  figure,  the 
agreement  is  excellent.  All  other  cases  examined  show  similar  results. 
Assuming  that  an  objective  of  the  estimation  procedure  is  to  develop 
the  best  representation  of  the  data  that  is  possible,  It  is  con¬ 
cluded  that  large  values  of  m(>  5)  are  not  necessarily  unreasonable. 


DENSITY 


FIGURE  1 .  A  COMPARISON  OF  THE  GENERALIZED  GAMMA 
DISTRIBUTION  WITH  EXPERIMENTAL  DATA 
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Also,  since  the  predicted  representations  always  represent  the  data 
very  well  for  both  large  and  small  m,  the  numerical  procedure  used 
to  solve  the  transcendental  equations  for  the  estimation  parameters 
is  judged  sufficiently  accurate. 

Further  investigations  revealed  that  in  some  of  the  cases 
the  predicted  design  extreme  values  were  unreasonably  high  as 
compared  with  the  most  significant  values.  In  order  to  determine 
more  precisely  a  reasonable  value  for  the  ratio  of  the  design  extreme 
value  to  the  most  significant  value  (XQ/Xs ) ,  an  analysis  of  over  100 
previously  processed  data  cases  were  examined.  The  results  indicated 
that  unrealistically  large  ratios  were  obtained  for  those  cases  for 
which  the  normalized  logarithmic  moment  (see  Equation  2.4  in  Reference 
2)  is  greater  than  -0.5.  Most  significant  values  were  realisitic  as 
compared  with  those  measured  in  trials.  Consequently,  in  the  current 
data  base,  cases  for  which  T  >  -0.5  should  not  be  processed  to  deter¬ 
mine  design  extreme  values  using  the  generalized  gamma  procedure. 

The  experimental  histogram  of  AR1  3E  30  has  a  substantial 
tall  which  may  contribute  to  the  large  design  extreme  value  predic¬ 
tion.  The  last  four  bins  (out  of  12)  each  have  only  one  observation. 

A  data  sensitivity  study  was  conducted  in  order  to  determine  how 
sensitive  the  design  values  are  to  the  tail  portion  of  the  measured 
data.  The  results  of  this  study  are  presented  in  the  next  section. 
SENSITIVITY  OF  RESULTS  ON  EXPERIMENTAL  DATA 

A  data  sensitivity  study  on  cases  AR1  3E  30,  S14  3H  10, 

SI 4  3Q  10,  and  S82  3Q  10  was  carried  out  in  order  to  make  an  assess- 
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merit  of  how  sensitive  predicted  results  are  to  the  experimental  data. 
Computations  were  carried  out  for  several  cases  for  which  four  obser¬ 
vations  from  a  rather  distinct  tail  were  dropped.  The  results  are 
presented  in  Table  1. 

The  results  indicate  that  elimination  of  a  small  percentage 
of  the  experimental  data  can  result  in  large  changes  in  the  design 
extreme  value,  particularly  when  T  >  -0.5  and  c  <  0.0.  Based  on 
these  results,  It  appears  that  a  substantially  large  number  of  obser¬ 
vations  in  the  tail  portion  of  the  experimental  data  is  desirable 
In  order  to  predict  results  that  are  insensitive  to  minor  errors  in 
the  observations. 

It  should  be  noted  that  the  data  sensitivity  study  as  conducted 
above  Is  rather  severe,  i.e.,  four  points  were  removed  in  a  systematic 
manner  from  the  tall  portion  of  the  measured  histogram.  The  proba¬ 
bility  that  four  points  from  the  tail  would  be  incorrectly  measured 
Is  rather  small.  One 'way  of  determining  how  long  an  experiment  should 
be  conducted  and  consequently  how  many  observations  should  be  made 
is  to  monitor  various  statistics  such  as  moments  during  any  particular 
experiment.  Once  these  statistical  measures  reach  steady  state,  the 
sample  can  be  judged  adequate.  Proceeding  in  this  manner  may  have 
been  impractical  during  a  single  prototype  run.  However,  combining 
several  runs  under  nearly  Identical  experimental  conditions  may  prove 
possible. 
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RESULTS  OF  SENSITIVITY 
STUDY  ON  FOUR  REPRESENTATIVE  CASES 
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EVALUATION  OF  THE  MAXIMUM  LIKELIHOOD  METHOD 


The  method  previously  employed  to  determine  estimation  param¬ 
eters  for  the  generalized  gamma  distribution  was  based  on  the  Stacy- 
Mlhram  method  (see  Reference  1  and  2).  One  of  the  advantages  of  the 
Stacy-Mihram  method  is  that  negative  values  of  the  parameter  c 
are  allowed.  However  based  on  the  current  analysis  of  over  100  cases, 
it  appears  that  cas-ss  for  which  c  is  negative  result  in  unrealistically 
large  values  of  the  ratio  of  the  design  extreme  value  to  the  significant 
value  regardless  of  tne  processing  technique  (Stacy-Mihram  method, 
grid  search,  nonlinear  least  squares  method).  Consequently,  development 
of  representations  of  the  experimental  via  the  generalized  gamma 
function  should  be  restricted  to  cases  for  which  the  parameter  c  >  0. 
More  detailed  analysis  (see  comments  at  the  beginning  of  this  section) 
has  revealed  that  processing  should  be  restricted  to  those  cases  for 
which  the  normalized  logarithmic  moment  is  less  than  -0.5*.  If  the 
above  restrictions  are  imposed,  the  preferred  processing  method  is 
the  maximum  likelihood  method  as  developed  in  a  previous  section  of 
this  report.  One  of  the  reasons  the  maximum  likelihood  method  is 
preferred  is  that  it  provides  unbiased  and  efficient  parameter  esti¬ 
mates  in  the  asymptotic  limit,  i.e.,  as  the  number  of  samples  becomes 


♦Examination  of  the  data  shows  that  It  cannot  be  adequately  represented 
by  any  of  the  standard  distributions. 
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large,  provided  specified  regularity  conditions  hold  (see  References 
6  and  7  for  details).  Also,  the  method  is  applied  to  the  actual  experi¬ 
mental  measurements,  i.e.,  grouping  the  data  into  histograms  is  not 
necessary. 

A  comparison  of  the  maximum  likelihood  results  with  experi¬ 
mental  data  is  presented  in  Figure  2.  As  demonstrated  by  the  figure, 
the  agreement  between  the  generalized  gamma  representation  and  the 
experimental  data  is  excellent.  The  ratio  of  design  extreme  value  to 
most  significant  value  is  3.65  which  is  less  than  the  value  obtained 
by  the  Stacy-Mihram  method  is  5.58. 

DIGITAL  PROGRAM  DESCRIPTION 

The  current  digital  program  is  based  on  a  previous  program 
that  implemented  the  Stacy  and  Mlhram  method.  The  following  modifi¬ 
cations  have  been  made: 

1.  The  maximum  likelihood  method  has  been  incorporated 
into  the  previous  program.  This  is  accomplished  by 
replacing  subroutine  GETCML  with  a  new  version.  The 
maximum  likelihood  method  should  be  used  only  when  the 
numbers  of  observations  is  greater  than  100,  or  the 
number  of  bins  is  greater  than  60.  The  Stacy-Mihram 
method  should  be  used  when  the  number  of  bins  is  less 
than  20  and  the  number  of  observations  is  greater  than  100. 
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2.  As  an  added  option,  data  can  be  read  In  the  form  of 
densities  or  counts  per  bin. 

3.  Printer  plots  of  Rayleigh,  generalized  gamma  function, 
and  experimental  data  are  generated. 

4.  Since  the  extreme  values  are  a  function  of  the  number  of 
observations  (time),  time  effect  plots  are  generated  along 
with  tabular  results. 

Three  proprietary  programs  supplied  by  DTNSRDC  are  used  in 
the  present  system.  These  are  PLOTPR,  MGAMMA,  and  MDGAM.  PLOTPR  is 
used  to  produce  printer  plots  of  data  written  on  unit  3,  MGAMMA  computes 
the  gamma  function,  and  MDGAM  computes  the  incomplete  gamma  function. 
Printer  plot  programs  and  programs  for  the  gamma  function  are  readily 
available  on  most  computer  systems.  Hence,  replacing  these  with  programs 
available  at  different  installations  should  present  no  major  difficulties. 
The  incomplete  gamma  function  can  be  obtained  from  the  gamma  function  via 
a  simple  integration  procedure.  However,  due  care  must  be  exercised  to 
assure  sufficient  numerical  accuracy. 

A  complete  description  of  input  data,  printed  output,  and 
subroutines  follows.  A  listing  of  the  program  is  presented  in  Appendix 
A.  An  example  of  print  out  is  presented  in  Appendix  B. 

INPUT  DATA 

A  description  of  input  data  follows.  Card  Images  of  a  typical 
run  stream  are  illustrated  in  Figure  3. 
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Card 

Variable 

Cols. 

Format 

Description 

1 

ic(i) 

1-2 

12 

option  to  convert  input 
data  from  frequencies  to 
densities. 

HNEW  =  x  DEL^ 

ON:  IC(1 )  =  1 

2 

IC  (2 ) 

3-4 

12 

option  to  convert  L  values 
from  time  to  sample/ run  time 

'"NEW  ■  (L0LD/Run  tlTO)  X  N 
ON:  I C ( 2 )  =  1 

3 

IC  ( 3 ) 

5-6 

12 

option  to  change  random 
variable  values  to  midpoint 
values 

ON:  I C ( 3 )  =  1 

1 

IC  ( 4 ) 

7-8 

12 

option  to  convert  input 
data  by  factor  CAL 

ON:  IC (4 )  =  1 

1 

IC  ( 5 ) 

9-10 

12 

option  to  plot  Histogram, 
Rayleigh  and  Gamma  with 
"PLOTPR" 

ON:  IC ( 5 )  =  1 

1 

IC  (6 ) 

11-12 

12 

option  to  plot  above  in 

MKS  units  (change  by 
factor  CONVER) 

ON:  I C (6 )  *  1 

1 

IC(7) 

13-14 

12 

option  to  plot  L  in  hours 
vs.  Xp,  Xq  (Time  effect 

plots) 

ON:  IC ( 7 )  *  1 

1 

IC  (8) 

15-16 

12 

option  to  plot  above  in 

MKS  units 

ON:  IC(8)  =  1 

1 

IC  (9) 

17-18 

12 

option  to  print  gamma 
functions  for  all  L  values 
OFF:  IC(9)  *  1 

ON:  otherwise 
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Card 

Variable 

Cols. 

Format 

Description 

1 

IC(10) 

19-20 

12 

option  to  print  above  in 

MKS  units 

ON:  IC(10)  =  1 

1 

IC(ll) 

21-22 

12 

scale  on  Time  Effect  Plots 

Each  dot  will  be  spaced  by 

10  **  (IC(ll)). 

-9  <  IC(ll)  <  99 

Default:  .05,  1-5  hrs. 

1 

IC( 1 2) 

23-24 

12 

IC(12)=0  calculates  design 
values  and  option  for  plots. 
IC(12)=1  no  design  values 
calculated  and  no  plots 
printed 

IC(12)=2  no  design  values 
calculated,  but  option  for 
plots. 

1 

IC(15) 

29-30 

12 

debug  printout  for  Stacy-Mihram 
method 

1 

IC ( 1 6) 

31-32 

12 

number  of  double  pages  on 
plot 

Default:  1  double  page 
minimum  is  1  double  page 
(2  pages) 

1 

IC(20) 

39-40 

12 

IC(20)=1 :  Use  maximum  likeli¬ 
hood  method  in  estimating 
parameters 

Otherwise:  Use  Stacy-Mihram 
method 

1 

ISTOPD 

41-50 

no 

Number  of  increments  in 
excess  of  200  to  be  used 
in  determining  various 
statistical  measures. 
Recommended  value  is  zero 
or  one. 

1 

IPLOT 

51-55 

15 

IF  IPLOT  f  0,  Cal  comp 

plots  will  be  generated 


I  RAY 
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1 


56-60 


15 


IF  IRAY  f  0,  Rayleigh 
distribution  will  be 
generated 


Card 

Variable 

Cols. 

Format 

1 

CONVER 

61-70 

FI  0.0 

1 

CAL 

71-80 

F10.0 

2 

I0P(1) 

1-2 

12 

2 

IOP ( 2 ) 

3-4 

12 

2 

IOP ( 3 ) 

5-6 

12 

2 

IOP (5) 

9-10 

12 

3 

N 

1-5 

15 

3 

K 

6-10 

15 

3 

ALP 

11-20 

F10.0 

3 

DEL 

21-30 

F10.0 

3 

CMNT(l-5) 

31  -80 

5A10 

4 

NL 

1-10 

no 

4 

ELL(1-5) 

11-60 

5F10.I 

X(l),  H(l)  1-80 


8F10.0 


Description 

Conversion  factor  to  go 
from  English  to  MKS  units 

Conversion  factor  of  input 
data  to  measurement  desired 

IOP ( 1 )=0:  c  =  2.0  as 
initial  guess  for  use  in 
Wegstein  iteration 
I0P(1)=1:  Initial  c  value 
read  in. 

I0P(2)=1:  Write  initial 
c  value 

I0P(3)=1 :  Write  c,  m,  X 
estimates  obtained  by 
maximum  likelihood  estimation 

I0P(5)=1:  Debug  printout 

Total  number  of  samples 

Number  of  divisions 

Value  of  \  -  F(x)  corres¬ 
ponding  to  the  desired 
extreme  value 

Division  size 

Description  of  the  set 
of  data 

Number  of  L  values  to  be 
considered.  Up  to  five 
can  be  used  (See  Equation 
2.15  in  Reference  2.) 

up  to  five  values  of  L 

K  pairs  of  ordinate  and 
experimental  density  values 
to  be  analyzed  4  pairs  per 
card. 
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Multiple  cases  are  processed  by  stacking  the  data  sets  with 
the  last  card  as  a  blank  to  serve  as  a  flag  indicating  the  last 
case.  However,  the  option  cards  (first  two  cards)  are  only  read 
once.  A  run  consisting  of  two  sets  of  data  is  illustrated  in  Figure  3. 


XXXXXXX, C  “  770  0  5  .  r  3  J  *1 ,  f  * . 

C«  4R0.E  ,  TTTT*  NNNNNNNNNN.CC  .K. 
4TT4CH.BIN45  «  BBBBBBBBBBB*  I  f>=  TTTT* 
4TT4CH,N$R0C. 

4TTACH.I9SL. 

UiSPT<t.IBsI*?t./N5P0C) 

B I  NIP  . 


Ml 

15 

0.01 

l’.O 

r?2  H  30 

5 

311  . 

P.  3  3. 

244'’  . 

499  8. 

ft  * 

.0075 

1«. 

.0147 

30. 

.021? 

*?. 

.0174 

54. 

.0103 

66. 

.0151 

7«. 

.ft  03« 

"O. 

.001  5 

102. 

.0005 

114. 

.000  = 

126. 

.0003 

1 7  P  . 

.0 

150. 

.0 

162. 

.0 

174. 

.1)00  4 

263 

14 

0.01 

3.0 

P02 

3r  10 

5 

262. 

54  6  . 

1  63S 

.  3?7=. 

645ft. 

1.5 

.01527 

4.5 

.  P’M  7 

7.  e 

.050  59 

10.* 

,04962 

13.5 

.0559(1 

16.5 

."3435 

19. e 

.04071 

2?. 5 

.  0 1  9ft  P 

25.5 

.01272 

20.5 

.00763 

31.5 

•003P? 

34.' 

.00127 

37.5 

.00127 

40.5 

.00254 

FIGURE  3.  CARO  IMAGES  FOR  A  TYPICAL  COMPUTER  RUN 
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PRINTED  OUTPUT 

The  printed  output  is  dependent  upon  the  options  that  are 
chosen.  A  general  description  of  the  printout  follows.  A  sample  run 
is  shown  in  Appendix  B. 

The  first  page  consists  of  the  user  options  chosen.  The 
options  are  only  printed  once  for  all  the  data  sets.  The  initial  guess 
for  c,  the  value  of  c  at  each  iteration  used  in  finding  c  and  the 
estimated  parameters  are  next,  if  maximum  likelihood  method  is  chosen. 

If  Stacy-Mihram  method  is  chosen  the  iterations  for  finding  m  are  given. 

The  second  page  identifies  the  data  set  being  processed.  The 
number  of  observations,  number  of  bins,  a,  and  bin  size  are  next 
printed  followed  by  the  bin  widths,  the  histogram  data.  Next  the  log¬ 
arithmic  moments  of  the  data  follows;  y,  s,  and  T.  The  last  items  in 
the  header  are  the  calculated  generalized  gamma  parameters  m,  c,  X  as 
computed  by  the  maximum  likelihood  method  or  Stacy-Mihram  method. 

Eight  columns  follow:  the  values  of  the  random  variable,  the 
value  of  the  random  variable  in  MKS  units  if  CONVER  is  given  and  Option 
6  of  the  IC  options  is  chosen,  the  corresponding  Rayleigh  predicted 
frequencies,  the  corresponding  measured  densities,  the  corresponding 
Rayleigh  predicted  densities  and  the  corresponding  generalized  gamma 
predicted  densities. 

The  values  of  L,  the  density  function,  the  cumulative  distri¬ 
butions,  and  corresponding  random  variables  are  printed  for  each  value 
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of  L.  These  occupy  eight  columns  of  output  arranged  as  follows: 
Column  1  -  running  index 
Column  2  -  value  of  the  random  variable. 

Column  3  -  f(x),  (Equation  1  in  Reference  1) 

Column  4  -  f(x),  (Equation  11  in  Reference  1) 

Column  5  -  f(x)^ 

Column  6  -  G(x),  (Equation  12  In  Reference  1) 

Column  7  -  h(x),  (Equation  14  in  Reference  1) 

Column  8  -  H(x),  (Equation  15  in  Reference  1) 

The  following  page  gives  the  Rayleigh  design  values  for  all  L  values. 

The  fourth  page  gives  the  generalized  gamma  design  values 

for  one  L  value.  These  include  x  .  xrt,  x_,  x„,  x„  ,  x_,  and  xn  .  To 

M  U  5  G  GA  u  DA 

the  right  of  these  values,  x_,  x_,  x„..,  x.,  and  xn.  as  computed  by 

U  u  uA  U  UA 

a  table  look-up  procedure  are  printed.  This  page  is  repeated  for 
each  l  value. 

After  these  pages  there  is  a  comparison  of  design  values  for 
all  L  values.  The  next  pages  consist  of  a  printer  plot  which  show 
the  measured  density,  the  density  as  predicted  by  the  generalized 
gamma  function,  and  the  density  predicted  by  the  Rayleigh  density 
function. 

The  following  pages  consist  of  a  time  effect  printer  plot 
which  give  L  In  hours  vs.  x^  and  xQ.  The  Cal  Comp  plots  generated 
replicate  the  printer  plots. 


PROGRAM  DESCRIPTIONS 
Main  Program 

The  main  program  serves  as  a  driver  for  the  entire  digital 
program.  In  addition,  many  of  the  required  computations  are  performed 
in  the  main  program.  A  description  of  program  flow  follows. 

All  data  required  for  the  program  is  read  from  unit  5.  A 
description  of  this  data  is  presented  earlier  in  connection  with  the 
input  data. 

After  reading  the  data,  the  moments  of  the  logarithms  of  the 
random  variables  are  computed  along  with  the  Rayleigh  parameter  and 
maximum  histogram  value  The  normalized  third  logarithmic  moment,  T, 
is  tested  against  the  value  -0.5.  If  T  >  -0.5,  the  case  is  not  processed. 
The  parameters  m,  c,  and  X  are  computed  by  a  call  to  subprogram  GETCML 
via  the  maximum  likelihood  or  Stacy-Mihram  methods. 

Once  values  for  m,  c,  and  X  are  determined,  the  various  statis¬ 
tical  measures  are  computed  for  up  to  five  values  of  L.  The  density 
functions  and  cumulative  distributions  are  printed  for  each  L  value 
considered,  followed  by  the  extreme,  most  probable,  and  significant 
values.  The  extreme,  most  probable,  and  significant  values  are  obtained 
by  the  Wegstein  iterative  method  (see  Subroutine  WEG). 

After  the  above  computations  are  performed,  Cal  Comp  and  printer 
plots  are  generated.  These  consist  of  comparisons  of  the  measured  den¬ 
sity  function  with  the  generalized  Gamma  density  function  and  with  the 
Rayleigh  density  function. 

Descriptions  of  major  FORTRAN  symbols  are  presented  in  Table  2. 
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TABLE  2 


MAJOR  FORTRAN  SYMBOLS  USED  IN  MAIN  PROGRAM 


FORTRAN 

Math 

Symbol 

Symbol 

Description 

N 

N 

total  number  of  samples 

K 

K 

number  of  divisions 

ALP 

a 

value  of  l-F(x)  corresponding  to 
the  desired  extreme  value 

DEL 

Ax 

bln  size 

NL 

- 

number  of  L  values 

ELL 

L 

array  of  L  values 

ISTOPD 

number  of  increments  in  excess  of 
200  to  be  used  in  determining 
extreme  values 

X 

X 

array  of  ordinate  values 

H 

f 

i 

experimental  density  values 

S 

S 

value  related  to  second  logarithmic 
moment 

YBAR 

y 

first  logarithmic  moment 

T 

T 

value  related  to  third  logarithmic 
moment 

R 

R 

Rayleigh  parameter 

C 

c 

parameter  c 

M 

m 

parameter  m,  real 

AMBDA 

X 

parameter  lambda 

XO 

X 

ordinate  value 

FOXO 

f(x;X»m,c) 

predicted  density  function 

BFOXO 

F(x) 

cumulative  distribution  function 

BFOXL 

F(x)L 

cumulative  distribution  raised  to 
the  Lth  power 

GOFX 

g(x) 

density  function  corresponding  to 
F(x)L 

HLO 

h(x) 

asymptotic  density  function 

HHO 

H(x) 

asymptotic  cumulative  distribution 

XSUBM* 

xy 

most  probable  value 

XSUBO* 

x 

0 

value  at  which  F(x)  »  2/3 

XSUBS* 

X 

s 

significant  value 

XSUBG* 

X 

g 

most  probable  exteme  value 

♦These  symbols  are  used  as  labels  in  the  printed  output  generated  by 
the  program. 


TABLE  2 

MAJOR  FORTRAN  SYMBOLS  USED  IN  MAIN  PROGRAM  (Cont.) 


FORTRAN 

Math 

Symbol 

Symbol 

Description 

XSUBGA* 

X9a 

most  probable  asymptotic  value 

XSUBD* 

xo 

design  extreme  value 

XSUBDA* 

xda 

extreme  asymptotic  value 

*  These  symbols  are  used  as  labels  in  the  printed  output  generated  by 
the  program. 
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Subroutine  GETCMl 


Purpose:  To  determine  the  parameters  of  a  generalized  gamma 
distribution  from  statistical  measures  of  input  data  by  either 
Stacy-Mihram  method  or  maximum  likelihood  method. 

Usage:  CALL  G£TCML(AT,U,E,S,YBAR, GAMMA, 00C,T) 

COMMON/PARAMS/./PRNT/, /HIST/, /INPUT/, /OPT/, /IOUNIT/ 


Subroutines 

and  Subprograms  Called: 

MLHCML 

Maximum  likelihood  estimation 

WEG 

Wegstein  iteration 

PSI 

Di gamma  function 

PSI1 

Trigairma  function 

DIST 

Evalutes  measured  and  observed 

frequencies  and  densities 

MGAMMA 

Gamma  function 

Description  of  Parameters:  See  Table  3. 

Remarks: 

(1)  If  Xp  does  not  converge,  a  check  value  is  sent 
to  calling  routine. 


♦ 


TABLE  3 

MAJOR  FORTRAN  SYMBOLS 


Name 

Location 

Description 

AMBDA 

/PARAMS/ 

A  parameter  of  output  distribution 

C 

/PARAMS/ 

c  parameter  of  output  distribution 

M 

/PARAMS/ 

m  parameter  of  output  distribution 

X(I) 

/INPUT/ 

random  variable  values 

ICK1 

/OPT/ 

Check  value.  If  ICK1=1,  c  root 
was  not  found  via  maximum 
likelihood  method 

rop(i) 

/OPT/ 

options 

YBAR 

argument 

logarithmic  mean  of  input  data 

S 

argument 

second  log  moment  of  input  data 

T 

argument 

third  log  moment  of  data 

CMNT(I) 

/INPUT/ 

header  label  for  GETCML  output 

FOBS,  FTH 

/PRNT/ 

computed  values  for  observed 
and  theoretical  frequencies 

DENO,  DENT 

/PRNT/ 

computed  values  for  observed 
and  theoretical  densities 

OOC 

argument 

reciprocal  of  c 

ICK 

argument 

check  value.  If  XD  does  not 

converge,  ICK=1  is  sent  back 
to  calling  routine. 

U 

argument 

trigamma  function  of  final 
m  value 

E 

argument 

diganma  function  of  final 
m  value 

TABLE  3 

MAJOR  FORTRAN  SYMBOLS  (Cont.) 


Name 

Location 

Description 

AT 

argument 

absolute  value  of  the  third 
log  moment  of  data 

H(I) 

/HIST/ 

densities  corresponding  to 
random  variable  values. 
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Subroutine  MLHCML 

Purpose:  To  estimate  parameters  X,  c,  and  m  by  the  maximum 
likelihood  method  for  the  generalized  gamma  function. 

Usage:  CALL  MLHCML 

COMMON/PARAMS/ ,/HIST/ , /INPUT/ ,/ IOUN IT/ ,/PRNT/ ,/OPT/ 


Subroutines  and  Subprograms  Called: 


UEG 

Wegstein 

iteration 

PSI 

Di gamma 

function 

F 

Function 

of  parameter  c 

G 

Function 

used  by  WEG  to  determine  c 

Description 

of  Parameters: 

Name 

Location 

Description 

AMBDA 

/PARAMS/ 

X  parameter  of  output 
distribution 

C 

/PARAMS/ 

c  parameter  of  output 
distribution 

M 

/PARAMS/ 

m  parameter  of  output 
distribution 

X(I) 

/INPUT/ 

random  variable  values 

0(1) 

/HIST/ 

densities  corresponding  to 
random  variable  values 

IOP(I) 

/OPT/ 

options 

ICK1 

/OPT/ 

check  value.  If  ICK1=1 , 
c  root  was  not  found 

Remarks: 

(1)  c  is  restricted  to  values  greater  than  zero. 


(2)  If  a  root  for  F(c)  Is  not  found,  a  check  value  is 
sent  to  the  calling  routine. 

(3)  An  initial  guess  for  c  for  use  in  WEG  may  be 
read  in.  Otherwise  a  default  value  of  2.0 


Subroutine  DIST 


Purpose:  To  calculate  frequency  and  density  distributions 
for  Rayleigh  and  generalized  gamma  distributions. 

Usage:  CALL  DIST  (all  parameters  and  arguments  are  passed 
in  COMMON  Blocks) 

Subroutines  and  Subprograms  called: 

MOGAM  evaluates  incomplete  Gamma  function 

Description  of  Variables: 


Variable 

Location 

Description 

X(I) 

/INPUT/ 

random  variables  values 

H  ( I ) 

/INPUT/ 

array  of  histogram  ordinates 
for  input  data 

K 

/HIST/ 

number  of  bins  in  input 
histogram 

N 

/INPUT/ 

number  of  observations 
in  input  histogram 

FTH(I), 

FOBS(I) 

/PRNT/ 

computed  values  for  the 
observed  frequency  distri¬ 
bution  and  theoretical 
frequency  distributions 

DENO(I) , 
DENTR(I) 

/PRNT/ 

computed  values  for  the 
observed  density  distri¬ 
bution  and  theoretical 
frequency  distributions 

Remarks:  The  input  histogram  need  not  be  evenly  spaced  in 
the  abscissa  x. 

Method:  Generalized  gamma  cumulative  distribution  is 
evaluated  by  MDGAM.  The  density  function  is  obtained  by 
numerically  differentiating  the  obtained  cumulative.  The 
Rayleigh  density  distribution  is  obtained  in  a  like  manner. 
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Subroutine  WEG 

Purpose:  To  determine  root  of  x  =  f(x)  by  Wegsteln  iteration. 
Usage:  CALL  WEG  (I,  X,  J,  N) 

Description  of  parameters: 

I  -  input,  iteration  count.  Initialization  action  is 
taken  when  1=1. 

X  -  input,  estimate  of  root  of  x  3  f(x).  On  output, 
refined  estimate  of  root. 

J  -  number  of  significant  digits  of  accurary  desired  in  x 
for  solution  —  used  only  in  initialization. 

N  -  output  -  completion  flag 

0  convergence  not  obtained 
1  convergence  is  within  given  accuracy 
Remarks:  The  function  whose  root  is  to  be  found  is  computed 
externally.  The  calling  sequence  is  as  follows. 

(1)  Initialize  by  issuing  CALL  WEG  (1,  XG,  J,  N)  where 
XG  is  the  initial  guess  for  the  root  and  J  is  the 
tolerance  described  above.  N  is  irrelevant. 

(2)  Set  up  a  loop  to  calculate  XN=F(X)  using  output 
X  from  WEG 

(3)  The  simplest  calling  sequence  would  then  be 

X*XG  Initialize  X  to  guess  value 

DO  1  1*1 ,ITERCT 

CALL  WEG ( I ,X,IT0L,N)  First  pass  will  inialize  WEG 
K  vN.EQ.l)  GO  TO  2 
X  -  F(X) 

1  CONTINUE 

C  MAXIMUM  NUMBER  OF  ITERATIONS  REACHED  WITHOUT  CONVERGENCE 
STOP 

2  CONTINUE 

C  ROOT  FOUND  WITHIN  TOLERANCE 
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Method :  The  Wegstein  method  refines  the  root  of  the  equation 


x  =  f(x)  by  calculating  an  improved  estimate  x^  which  is 
intersection  point  of  the  line  y  =  x  and  the  secant  line  of 
f(x)  based  on  the  previous  two  evaluations  of  f(x).  This 
method  requires  only  one  function  evaluation  per  iteration. 
The  equation  for  the  intersection  point  is  given  by 

t  =  FU^  *  *x-'  ~  FC**-~>  » 

Completion  code  N  is  set  1  when 

iXXH  -*,1  <  1  I01*  XxJ, 

i.e.,  when  the  change  in  X  between  iterations  is  less  than 
one  part  in  10^. 

Function  PSI1 

Purpose:  To  evaluate  the  trigamma  function,  f(x) 

Usage:  PSIl(x) 

Description  of  Parameters: 

PSIl(x)  -  value  of  the  triganria  function  at  argument 
X  (Output) 

X  -  argument  of  the  function  (Input) 

Remarks: 

(1)  Argument  x  must  be  greater  than  zero. 

(2)  The  trigamma  function  is  the  second  derivative  of 
the  natural  logarithm  of  the  gamma  function. 


the 
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Method:  For  arguments  greater  than  13  an  asymptotic  expansion. 


is  used.  The  truncation  error  associated  with  the  above  expan¬ 
sion  is  less  than  1  x  10"^  for  x  >  13.  For  arguments  between 
0  and  13,  the  recursion  relation 

^  Cx)  =  ^  l  1  - — 1 - rx 

i*i  Cx  -ri  -O 

is  used,  where  m  is  chosen  as  the  smallest  integer  for  which 
x  +  m  >.  13.  'fl(m  +  x)  is  evaluated  by  the  previous  formula. 

The  formulas  are  obtained  from  Reference  8. 

Function  PSI2 

Purpose:  To  evaluate  the  tetragamma  function,  '{'"(x) 

Usage:  PSI2(x) 

Description  of  Parameters: 

PSI2 (x)  -  value  of  tetragamma  function  at  argument  X 
(Output) . 

X  -  argument  of  the  function  (Input). 

Remarks: 

(1)  argument  must  be  greater  than  zero. 

(2)  the  tetragamma  function  is  the  third  derivative  of 
the  natural  logarithm  of  the  gamma  function. 


Method:  For  arguments  greater  than  13,  the  asymptotic 
expansion, 


\ 

'of 


is  used.  The  truncation  error  associated  with  the  above 
expansion  is  less  than  1  x  10"^  for  x  >  13. 

For  arguments  between  zero  and  13,  the  recursion  relation 

r  go  =  -  i  , 

is  used,  where  m  is  the  smallest  integer  which  satisfies 
x  +  m  >_  13,  and  H'"(x  +  m)  is  evaluated  by  the  previous  formula. 
The  above  formulas  are  obtained  from  Reference  8. 

Function  PSI 

Purpose:  To  evaluate  the  digamma  function,  ¥(x). 

Usage:  PSI(x) 

Description  of  Parameters: 

PSI(x)  -  value  of  digamma  function  at  argument  X  (Output) 

X  -  argument  of  function  (Input) 

Remarks: 

(1)  Argument  X  must  be  greater  than  zero. 

(2)  The  digamma  function  is  the  derivative  of  the 
natural  logarithm  of  the  gamma  function. 
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Method :  For  arguments  greater  than  13  an  asymptotic  expansion 
is  used: 


u* 


A*  V  *  lAP**'”' 


The  truncation  error  associated  with  the  above  expansion  is 
than  1  x  10"! 0  for  x  >  13.  For  arguments  less  than  13,  the 
recursion  relation 

<PCx)  -  'pUt™)  -  4  -1 - 

(el  ^  ~  l 

is  used,  where  m  is  the  smallest  integer  which  satisfies 
x  +  m  >_  13,  and  y(x  +  m)  Is  evaluated  by  the  previous  formula. 
The  above  formulas  are  obtained  from  Reference  8. 

Subroutine  MGAMMA 

Purpose:  To  evaluate  Gamma  function 
Usage:  CALL  MGAMMA  (X,  GAMMA,  IER) 

Description  of  Parameters: 

X  -  Input,  argument  of  Gamma  function 
GAMMA  -  output  value  of  Gamma  function  at  X 
IER  -  error  code 

Remarks:  This  is  a  library  routine  supplied  by  DTNSRDC. 


User  must  attach  library  (IMSL). 
Error  code  meaning  is  unknown. 


Subroutine  MDGAM 

Purpose:  To  evaluate  incomplete  gamma  function 
Usage:  Call  MDGAM (T,  X,  PR,  IE) 

Description  of  Parameters: 

T  input,  limit  of  integral 
X  input,  exponent  in  integral 
PR  output,  value  of  incomplete  gamma  function 
IE  output,  error  code 
Remarks: 

(1)  This  is  a  library  provided  by  DTNSRDC 

(2)  User  must  attach  IMSL  library 

(3)  Meaning  of  error  code  unknown 

(4)  Function  evaluated  is  given  by, 

PCx/t)*  TTx)  f  o’t"<-'u<lu 

O 

Subroutine  HISTO 

Purpose:  To  set  up  labels  and  plot  the  experimentally  measured 
density  function. 

Usage:  Call  HISTO 

Subprograms  and  function  called:  SCALE,  AXIS,  SYMBOL,  LINE 
Remarks :  This  program  takes  the  measured  density,  which  is 
passed  in  COMMON  and  uses  these  values  to  generate  data  for 
a  Cal  Comp  plot. 
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Function 


Purpose:  To  be  used  by  WEG  in  evaluating  the  c  parameter 
by  the  maximum  likelihood  method. 

Usage:  G(c) 

Subroutines  and  Subprograms  Called:  PS I 
Description  of  Parameters: 

G(c)  -  value  of  function  (Output) 
c  -  argument  of  the  function  (Input) 

Remarks :  c  must  be  greater  than  zero. 

Method :  To  implement  the  Wegstein  technique  the  equation 
f(x)  *  0  must  be  put  into  the  form  x  =  g(x).  The  equation 
used  is  given  by 


c.  — 


- ■  [  Uc  f  \a  (.1  *i‘  °0 

-v  w  -  u(|o;) 

u  i.  *• o  •  — 1 -  >■  ' 


+  V 


r 


C 


1 

4  ~ 


i*' 

T7~ 


0i  Is  the  density 


Is  the  value  of  random  variable 


Function  F 


Purpose:  To  evaluate  the  function  used  to  find  the  c  parameter 
by  the  maximum  likelihood  method. 

Usage:  F(c) 

Subroutines  and  Subprograms  Called:  PS  I 
Description  of  Parameters: 

F(c)  -  value  of  the  function  used  to  find  the 
parameter  c  (Output) 

c  -  argument  of  the  function  (Input) 

Remarks:  c  must  be  greater  than  zero 
Method:  The  equation  used  is  given  by 


F(c)  =  Inc  +  U  (£  *oc) 

v  fi 

5  V  o.- 


W; 


£  o;  J 


—  lv\  k  o 


,  >< 

* 


—  c.  o;  ^ 


=  0 


< 

f.0; 

<U 


0^  is  the  density 

is  the  value  of  the  random  variable 


CONCLUSIONS 


An  extensive  study  of  an  existing  digital  computer  program 
for  determining  parameters  required  to  represent  SES  1 00B  trials  data 
by  a  generalized  gamma  function  was  carried  out.  Difficulties  encountered 
in  earlier  studies  in  which  approximately  5  percent  of  the  cases  analyzed 
produced  unrealistic  values  for  the  distribution  parameters  have  been 
resolved.  In  order  to  resolve  this  problem,  data  sensitivity  studies 
were  carried  out  and  three  other  techniques  besides  the  Stacy-Mihram 
method  were  investigated.  Based  on  results  of  the  study,  criteria  are 
developed  for  determining  cases  that  should  not  be  processed.  A  summary 
of  major  conclusions  follows. 

a.  The  Stacy-Mihram,  maximum  likelihood,  grid  search,  and 
nonlinear  least  squares  methods  all  yield  estimated  proba¬ 
bility  distribution  functions  which  represent  very  well 
the  experimentally  measured  histograms. 

b.  Design  extreme  values  are  more  sensitive  to  the  experi¬ 
mental  measurements  than  to  the  numerical  procedure  used 
in  determining  the  estimation  parameters. 

c.  In  some  cases  the  predicted  extreme  design  values  are 
unrealistically  large  whereas  the  significant  values  are 
realistic.  Based  on  an  analysis  of  over  100  cases,  it 
Is  concluded  that  cases  having  a  nondimenslonal  third 
logarithmic  moment  greater  than  -0.5  produce  unrealistic 
extreme  design  values.  Results  of  a  data  sensitivity 
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study  have  revealed  that  in  all  these  cases  the  data 
may  not  be  considered  as  a  sample  taken  from  a  steady- 
state  random  phenomenon,  and  hence  it  is  recommended 
that  statistical  analyses  for  these  cases  not  be  carried 
out. 

d.  The  Stacy-Mihram  method  should  be  used  when  the  number 
of  bins  (number  of  divisions  used  in  the  analyses)  is 

small  (<  20)  and  when  the  number  of  samples  is  large  (>  100). 

e.  The  maximum  likelihood  method  should  be  used  when  the 
number  of  bins  is  large  (>  60)  or  when  the  sample  size  Is 
large  (>  100). 
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APPENDIX  A 

COMPUTER  PROGRAM  LISTING 


PRO'-RAK  ’?H  (T'lPiJT  ,0'jTPLT,TAP£5=Tl'!Pl!T,TS«'  r£  =  0'JT°UT.TARP7.TAPr3) 
£s  L  waxx 
’£AL  * 

•'THESS  !•"■  m#«00)«2»*90».F,l,»e90)  •  “»: C  rCfl  J  .  1 7 « 0  ( 4  ) 

3IS£>SIO':  X*;''S*5'jU  (500)  .FNF.jp  oo  ).  ’v(c  00) 

T'E'tSn*’  Y0vS(?"  j  >  .  pry  "'■» (  2TT  ),  -fo  *  *00  >  .8F0L(2'j  3)  . 

*  Snr  X-''(20''>»FLO“S(2oa>,H,-‘0'*,'*',OG) 

?HEMSI0'  PDX(7).ELT  *7)  .FSC(7>  »PPX««(  R)  .F'-Xm  .B,'X"(7) 

CC5Y0N/1’  °ut/x»100  )  ,rLL  (7>,'I,«LC  ,7(1  00  )  ,  ISTCP')  ,  I  PLOT  ,  I  =  »  v,  jrr?n  j 
COMMON  X0cT/I0P(23>  ,ICKt 

C0u“0»J/PL0/r.'PUN«fliPLC  T.rTP(<),lTY*5>.7TX(x>. 

V"J'*°  *  3  SXSCS  ,SCA(  10  >  «F  =  '•"’(  10) 

CO**  »  ON  .'PPMT  xrrp^<l  0  3  )  »F  Tu<  IOC)  .C5'lr  (  1  00)  .05r.T(13  0) 
i,FTH;»{100).Ur'iT',<iaO).C'",yP':,C4L 
COM'*-^ /»  IS  T2 /FIRST  X.OELX.FIROT  Y,  RELY 

'O'1 '•ON  IST/PIf.  LO,K»P  IV.’U  V)  ).«{  10'’  )  ,  Y  “  A  x  ,  XLF  *!.  YLF*' ,  Cw  T  («  ) 

COX'*' 'I /°4Ri  A  MB  04  iC(' 

701 » ON  /’P'JNIT/IDP^IFO.IVPT.IPCH 
OATS  JOB  ,1  =  "  . IW°T,  IPCH/t.E.6  .£  / 

:.?7A  XLr‘ .YLEN'l  j.«-!./ 
tTu=0./t. 

'•ALL  °1‘  75(3.3.7) 

^  Jrj.j  nor  15*1^ 

’£40*5.1’)  (  ’C*  I),I=l,2:),ISTOPn*IPl.CT,lP4Y,nnVFB.C4L 
«rITr(C«1413)  ICL?T  *  I  OS  v 

3413  ! C  OPTIONS  usrr  IV  TOTE  9NN  eft  tdlOTp.I’.xv. 

V=Hl04Yr,T=i 

00  1*  Irl.oo 

14  IF*  K(T)  .CT.O)  JOiTCf,;,  ]»!*»  I.ICH) 

1414  =}94ST(li  ,-lu  PPT 1  r  ,:?,1U-,I2» 

J°ITE(F , '313) 

2017  P09*'4T  (lF0,3?MJ(>o  QP  T  IC  *  c  US5C  T  *'  this  P'jN  SrT) 

’040*5,114)  ( I pp  *  I )  .  I  si »  *  0  > 

115  P07VST  f2rI?) 
n0  <379  I  p1  * 70 

Ic<  IPP(I  ).ft.o>  J9ITE(t,l4i4)  (I.nc*!)) 

970  COMTI'JIE 

JF(  CO*IVro.TQ.:  .  )  CO»)V5°sl.O 
15  FPR',fT<2''!?«71C»2I5»?Fl,'.&> 

<* 

09  CONTINUE 
»\« 1 )=T. 

C*K  1)  =  0. 

TP=1 

1 )=C. 

<LL  =  XLF‘)«? 

C 

C  •  •••*  R  E  AO  INPUT  SNC  PRINT  ..... 

C 

’54  3  (5.100)  ".R.ALPtDEL,CM'' T 
Tc»x.€0.0)  0 0  TO  10303 
no  rCR 'AT(.?I5.2c10.0,5A10> 
c 

C  ’51 J  L  VALUES 

’£40(5. 1"2>  4L»(5LL(  I).  I  =  1.'IL) 

13?  FORMAT  (I10.7F1C.0) 

75*0*5.1 ?1) (XCI) ,H( I ).I Pl.X) 


on  ji-g 


131  FORMAT*  =  F10.;  1 
230  CCR  “AT*  lx.  t'10  .*  1 
=  oa  =  l.-4L° 

C 

C  !r  UNEQUAL  I’.’TPPVALS  OrLr  1 
C  OEL  RE43  I*.  t  A RMY 

C  Iff  ST  4  ST  TM  F  0T!T  0.5  -  SIGN  JILL  PcrrrEr' 

C  THE  STARTING  VALUE  GIVE"  AS  CEL. 

OuX=DEL 

IF*OEL.rT.g.g)  =0  TO  1=0 

ST=-DEL 

0*1  l=X  (1  ) 

0SUM=0  <1  ) 

«(1):S’ 

CPXsC*  11 
0  =  1  =  6  T=2.* 

2*1 )  =x « i > 

!=< =< I 1.2E.0 *1-1 ) >  ?"Vrn<T) 

V  C  )  =X  <T  -l>»f'<  I  -  1  > 

? IMW ( 1 )=C< I) 

OSU-,=OSU,'*C<1> 

I  =  6  SCNTl’JL' 

Tc<  !C(  2)  .E'.Ol  TCLG  =  1 
TF< IFL=.E2.1>  G=  TO  1000’ 

'0  TO  1=0 
iso  =0  1==  T=1.X 
1=6  =  *1  >  rO  EL 
l  =  o  2"'1T1‘ILE 

C  1*>TT3N  I!  CONVERT  INPUT  T=  CE’ISITY  IF  •ECESSAcy 

:f(  icd)  .ne.h  ==  t=  = 

00  II  1=1,' 

II  J*I)  =  M<I)  /  *=LOAT*Nl«=tI>) 
r 

c  optic’;  2:  ’oorcr  l  values  if  necessary 
o  IF*  !2<?>  .ME  »  1 1  "•=  T=  I 

”Ti“E  =  ell  n i 
00  12  !  =  1*>'L 
r  LT  *  l )  =  ELL*  I) 

12  ELL  *  I )  =<ELL'T)/STI''E)»’ 

C  2°TI3N  I  :  4Q  JUST  I*'PUT  X*I)  TO  «irP=INT  ,0=  INTERVALS  Ir 
5  IP*  IC*?>.NE.l)  SC  TO  * 

20  n  1=  1,  x 

i *  mi  =  mi  *  oiri/i.s 

C  0POI0N  <t ;  CONVERT  INPUT  TC  "rs  SURE ‘'ENT  ’'ESIREE 
A  IF* I  C  <  4  >  ,Nr  •  l  )  SO  TO  = 

00  7  I  =  l.K 
U*I 1  =  “f II  *  CAL 
OtIl=D*I)*CAL 

x « i >  =  x *n  •  cal 

CONTINUE 

IF  L  VALLES  ARE  OPAD  I*1  AS  SABLES  =E»  TIwr  I‘TE»v»L 
COv=UTE  TI VE  VALUES  FOR  PL^T 
IF*IC*2).E0.1)  =0  TC  a 
RTI'»E=FLL*1)/,I 
20  ST  I  i  l.AL 
ST  ELTTT1  =  ELLft)  •  PT'E.'” 


nno  o  O  f) 


& 


c 


c 


continie 
a=: . 
yka  x=o  . 

CH=  1  . 

*>1  =  0. 
?2  =  C  . 
S3  =  3  . 


COupUTE  P,  »T 


00  1  I=1.K 
'1  =  H<t)*',<T) 

IF(T!.*T.Y>'A'')Y*‘AX  =  TI 
p  =P  *T1  »X  f  T  )  •*?. 

YJrALOGfYd)  ) 

<;h=<;h.»ti 

vr?-»i  *Y  t 
YU  =  71  •  *  !2 
Cl=? I ♦ T1  »YT 
S2=S2-.T1 *“!  2 
S3=S  **T!  »ytt 
continue 

.....  COv=UTF  YPS°.  S.  T  ..... 


rv=>j 
Y04°=S1 
TlsCM/  <E 1 .  » 

T?sTl*E»V(r»|.2.) 

T;=s? 

Turvgjo.yeio 
t=-T4.  '*w 

r=Tl»<T?-2..T4.T“) 

T:  T2.fS3-I.»YEi°.T3.YPe».f3..T»-T<:)  )/C.*l.I 
IFf  T. ST.  f-.G)>  UPITEf6.123Y1 

1237  FOP',«T(l>'0,Pf.MT  IS  .“EATER  THAN  -.=  '}'•  '*•  T"  NC  Y  T  ?»TA  rET) 

Ic  ^0  T "  P’ 

ATsASS  <T) 

.....  CCI"PUT"  M.r»LA'*oOA,  AMO  GA*“A  ..... 

I  -n 

CALL  GETCwLIAT.!l.E*S.YB?P.GAM*,A.r'nC.T.IC*r) 

IFCICK.EG.1)  GO  TC  i" 

IFC ICKl.fO. 1 )  GO  TO  09 
V9tTF(<,P01«J 

p  0 1  A  “39  M  ST  11  hj  ,  2?UP  *TLEI  C-H  0ESI6N  VALUE1-) 

00  “012  0=1. "I 

X®0=SQ9T  (P.ALOGf  rLLt.J)/  >-LP)> 

XP3MS=X9r.C0AVr9 
JPTTEI6.P013)  FUUd)  .X9 O.XPDvS 
»n?  C  0*1  TI»I  OF 

8313  FORMAT'!*  ,3PL=  «F1 0  .4 ,  =  X «RHXO SUEOs  .  FI  0  .  A  ,9  Y  ,  1  ?“XBSU0''  I  W“'S  )  =  , 
TF13 .A) 

TF<  ICO  I?  )  .FO.  1  >  BO  T0  1  ? C  9* 

c  .....  COMPUTE  C»v.LA«=0A  »ROO"CTS  Fno  LATER  USr  ..... 


V. 


T 1  =  C*" 

Tlv  =  Tt-1. 

T2  = 

T  'i  -  AEEf  T7) 

C>'RPC  =  Tl,‘«OOC 
r 

C  .....  10000  LOPO - OVER  Tur  L  VALUES  ..... 

TO  10000  IL=1.NL 

C  ..  INITIAL  GUESSES  FOR  X  AND  X  .. 

C  r.  0 

TXD  =  1000. 

E/D  =  0. 

EX3  =  0. 

TX3  =  C. 

TXO  =10  *0  . 

E  v0  -  G  . 

TX3  A  =■)  . 
rx3A=0. 
txo a  =i  :ao. 

'XQAIO  . 

EL  =  E  LL  < I L  > 

,°Irr  <0.10001)  EL 
11001  c  OR  *  AT  <lwl  *3uL=  ,  F  7  •  1  '  ) 

IF  <EL  .IT.  J10)  GO  TP  °001 
r'-'L  =  E»C<-CL) 

=  0  To  a  «  n  j 
R001  E  VL  =  •:. 
r002  C0'lTT»Ujr 

IF<IC(E)  .EC.C)  WPITEIfi. 004) 

TOO  =01  xi  jt  •  xy,.!»t=x,.x*»l  *,12X»«F  t  X)  X)L*. 

1  .1<xw, IPX, *LH<X)*.1'X. *«</)♦) 

770  FORMAT  <//»  EXTRAPOLATED  VALUER  »//> 

I  STOP  =  200*!ST?PO 

C  »•  SET  UP  “ISTPGRt'"  P»DA“FTE»S  ». 

I  =  7 
J=0 

a.INLP  =  x<l)-D<l)/2. 

-<AX  X  =2  .n*'X»xLE’l 
X0=  X<1)-D<1) 

<oo=xo«n<i )  /  2. 

IF<  X  0  D  .LE.  .31)  X0=0. 

711  TRESrlX? 

RES = <FL1AT< I )/?.)- I»FS 

IF<  C»ES.E'J.O)  <J  .LT  .x) )  J=J*1 

SELL  =0  <J>/2.0 

1FI'ir  =  C<v')/10.0 

XOsXCtOR  ’NC 

r  =  i  ♦  i 

IF  <1  .RT.  TSTCF  )  '-0  TO  1 1002 

XTO  =  /P 


C  »»  11010  L30°  .» 

00  11000  IF=1,E 
C 

C  .«♦»«  CONFUTE  S*  ALL  C'X> »  P«0oAf!LITY  0E'lsrTT  =*! 

c 


xo  =  XT7VJFIRCMIF-1  ) 


T4=XD»*TlM 

T45=(4V9DA.X'>»»C 

Tj  =  frvn<-T4e> 

-0X0  =  T=.T4.T5 
C 

Ic( IL.C-T.l)  rr  yo  10002 
IF< X?.3T.'*ivy 1  to  TO  10CC2 
I  p( IP. £0. 5C0>  GO  TC  10002 
C 

IP=IP*1 
=  *t  1  =  1  sFOXO 
»*.<  IP1  s»? 
xTlsXO/N 
XT2SX0  .x  T1 

=  NX  =  2.  .X  T1 »  =  xo<  -XT2  > 

:*'l  f  p  j  r°  X 

I(r<  «‘tx  .OT.  VMS  X)  ywAXs= ‘ X 

r  cf  e-x  r.r,T.  Y«iX  IXMAXsFO  *0 
103  02  CON  T IN  Lr 

C  .....  COv«»UTE  C*=ITSL  f  (»  1.  ==C=t°ILITX  “'I5TNINL1  TI 

r- 

CALL  MCOA'-f  T4<i,'*,P0.  Trl 

Ic  (IE. E-.12C.fP.IE.E0.no)  UNITE  («r,sm  IC 

531  r  ON  MAT  S5qH  F  =  NCN  IN  «»*•*(>  0 1  STN  I  °  IJTJ  ON  CL’MTTION  CO-ouTSTI 
.Ns  .181 

IP  fC  .LT.  0.1  ON  s  l.-BO 
=  P?XC  s  =  ° 

c 

TTT  s  tL*U1'(,!C,1l(5) 

Ir  f TTT.LT.-?03.)  •"  TC  =  "03 

?F3  XL  s  c'P'»yc»»rL 
CO  TC  “304 
='03  cFOXL  s  '. 

N"3»  '‘ONTIMLE 

C  .....  CO^NUTE  'XjlL  C-tX),  C**»LL  H«Y1.  C»°ITFL  “(XI 

C 

!F<  IC< 12  1.E3.21  r0  T 5  11000 
OOFX  s  EL*(?FCXl/BFOX3> .FOXO 
F X  s  EL* (ON-1 . ) 

I c  (EX  .LT. -300.1  GO  T"  =010 
EX  s  EX°  <FX  1 

UH()  s  EX-F'*L 
NO  To  8311 
=310  HLO  s  0. 

HHO  =-E  »L 
=011  C3NTINLE 

IF  (TF  .GT.  1)  GO  TO  11C01 
XOXS  (I  1  s  xo  •  CCNVrD 
FOXOMSU)  =  FOX?  /  CONVrN 
9F0  (  II  s  FFCXO 
NFOL(I)  S  =  F  n  XL 
50FX“(1>  s  GrFY  /  CO'IVE® 

HLOPS(I)  s  HL0/C2NXE® 

HHOrS(I)  s  HHO 


ICON  =  I 

IFI  ICC  *>  .£".1  >  TO  31  001 

UPlTf  •‘-.,710  T.XO.FOXo  .£Fr'XO,=FOvL,''CFir,HLC*HHO 
715  F09'*AT<lX.I,=  .7':i9.6> 

tiom  continue 

c 

c  .....  iniTI4LI7»T70'! 

c 

TX=AFS  IB^Xl-O^A) 

TF  <  TX-TxO>  o00C*900t,°001 
9"00  TxO  =  TX 

5 X 0  z  X7 

•=ooi  r.oNrr‘i(jr 

!e  fTXG.3E.~:FX>  GO  TO  a0?2 

txosoofx 

~v3  =  »0 

90  02  5CN71NLE 

Tx=ABS  '"OX^-TTh) 

TF  fTx-TXO)  95  3  3  *99  5  4*  "  3  4 

95  33  7<q  =  TX 
E*0  =  XO 
T"  34  OONTINLF 

!F  (TXOS  .^F.  ML?)  "0  72  °2  59 

’■X3S  =  uLC 

FXS4  =  XI 
90  09  CON  T  TN  IF 

TX  s  4eS<H“2-0«t> 

TF  (TX-TX04)  95:6*9*«CT.^0CT 

9500  Tx  0  A  r  *x 

;x3A  =  »  n 
95  37  CONTINl9 
11  C  00  C0N7 r *•  L “ 

X5=XT0*nrLL 

TF  f?F0X0.LE..9°9999 >  o"  T0  701 
00  TO  11053 

1  10  02  V 9 1 T "  <6*11504>  ?ST",D 

11004  FO»MAT  fl***90.  OF  IFCOIOS  5-T.1X.T*) 

11003  CCNTI*HE 

1FC ICC  10  >.NF.1>  GO  TO  17 
W«I TFf  £,7111 

711  FCT" ST(/ <lvtT40 .13H9 .x.5.  5YSTEV //) 

J9ITFtf.30*l 
00  1G  «.=  l.ICCI 

,9I  TF1S,71C  ,  J.XuMsr  J»  .  90X5l*3<  J)  ,9F0  1 J)  *cFnLC  J»  ,OPFv-(  J> 

*  “LC'S  *d)f  Hfrvisc  J) 

19  COXTINLF 

17  CON  T  TINLE 

C 

IF«IC(221.E0.2>  r. 0  70  °70 
C 
C 

C  .....  SOLVE  c09  X  •*»•* 

c  ** 

c 

79=1  «-COC>  «.OOC 
X  SU  B  v  =  T9'A»F0A 
X5'J9',S  =  XS1;0"  »  CO'. VE9 


W9IT-<S,2C7)  »SUPH,*Slfn,*S 


c 


73 


C 

C 

2os 

20F 


33:0 


c 

c 

c 

r 

C 

r 


C 

c 


3? 


702 


93? 

*3 


•  solve:  * c r  v  usi*.,r  yrssTri'  ite9ative  *ET“"r 

n 


T3=0“A»*a./EL) 

TT=AL°  >F  L 
U  =  rXC 
•IE=0 

CALL  VFGri.U.4.0) 

..  s  T  A  9  r  t’-rp  irn“  LOC9  «. 

T95  =  <A««rA.U>.»C 
CALL  «r'-i"MTA5.».PS,  TE7 
IF  (C  .LT.  O.J  99  =  l.-99 
J  s  J*99  -T-» 

ms*J*TT  *'L09CP9) 

CALL  WEC<f2«U,0.»C> 

JFITE  (A.?’?)  U.Pf*,T-».TT 

c09«*AT  (lX**Lr  *.023  .R.cv,»««  ,G  2". 9  ,k».«T7.,  520  .«  .5X..7T* 
IF<U  .17.  0.1  Us.000  1 

ne  s  *it  «  i 

IF  <»£  .LT.  =0)  SO  70  ’CSC 
J9ITE  <*.?'■$) 

rC9  x  A  T  A  7  ak  *'9  CO'.'VE’OC'C9  99  TA  r**rO  F  ?P  x9l«SJ 

•j  =  0. 

SO  7377 
CON  T  J  »it£ 

IFXNC.EO.OT  GO  79  3 3 

..  END  IT  rF  ATT  ON  cn9  *  •• 

VTJPO  =  U 

«.*»»  EVALUATE  C-A*»A  DIST9I9UTIi^  FU'CriO*'  ..... 
*Tr  (A'*=OA.'J).»C 

call  ocr-A^c  tt.^.pf.te;) 

.’ 9 1  7 F  '5,707?  Po 
9  39  M  A  T  »•  09  .,  "20  •  9  > 

.....  SOLVE  r0P  <  US  T*!S  WFOSTFIN  ITF9jrivr  «rT((?r 


CA  r  C*A-eOA 

OVC^U  ,-c»v 

C«3  =  C*»-7. 

EL*  =  EL-1. 

U0  =  £»0 
'J  =  uo 
NE  =  0 

CALL  iircci, U.A.O) 

*.  STA9J  ’TFB  A  TION  LOO 9  ♦* 

•J=0 

•U  -  A“onA«U 
AUC  =  AU**C 

Trl  AUC  .LT.  20.1  SO  7"  S001 

'J=l  .F-l" 


'0  TC  6’jZ* 

4001  COM  T  IMLf 

CX°T  =  F  vo  (  -4 UC  ) 

F=£Lv«(Ta«FyoT*'l»*Tl“)»*'> 

CALL  M3GAfMAUC,*',PR.  !EI 
IF  (C  ,LT.  0.1  PR  a  I.-FP 
]=P  /  <  <  fs  *U*»CU?«FXPT»  »  ONCw*C*AUC  >  >*TC  J 
5000  CONTINUE 

CALL  VFr-.<2,u,0«NCI 

:  W ® I  T r  <5, 7C3J  U.PR»F.C'<7,EXPT*C'MC',tAUC«Al'»T? 

703  FgRMsi  (,  orpije 

TF  (I!  /.F,  0.1  r-0  Tfi  6  0  0  3 

,JO  =  U'!*C<  11  O 
■I  =  (10 
OO  TO  »3= 

60CT  CONTIMt*- 

*;  r  —  »ip  *  i 

IF  <“€  .LT.  COT  SO  ”  *000 
JOITF  <6  .20  c  ) 

«*»a  =or«jt<x«w  NC  CCMVERSENCr  09TAINED  F nD  **U?SJ 
•  j  r  0 
*,n  TO  A- 
4000  COMTifa. 

IF<  '•C.EQ.G  )  CO  T9  43 

r  ..  fmq  ITFOATIQM  Fno  y  >. 


4  0 
437 


/$J  90  U 
CONTINUE 


SOLVE  F!'R  V  US  I  MO  VEGSTEI  *■  ITEPATIVE  "ETHOS 


C  'JO  -  X  ( f  ) 

"KsX  /3 
U0=  X  C  VK  ) 

t!  =  U3 
"E  =  0 

FALL  WES  <1.L,*»0> 

C  ..  ST5°T  ITFPSTJ0*|  LOOp 

5007  »U  =  AXaCA«U 
ADC  =  tU««C 

call  '«cga'« r alc .“»pr«  ie> 

IF  (C  ,t-T.  0.)  PR  =  1.-"’ 

1 1  -  U+TTL'— PR 
CALL  WEG  f2»t)40«'!C) 

C  VPITE<6xT03>  U«FP 

IF  (U.GE.O)  GO  TO  603* 

u=ja*ou ) 

C  n  =  uo 

D=l  .E-13 
600*  CONTINUE 

>’£  s  ME  «  1 

IF  INE.LT. SO)  GO  TO  60  99 
WRITE  (6.20°Q) 

2000  f  or  m  a  f  f  x  am  "0  CONVEROFxiff  OBTAINED  FOR  X  T’JB  0 ) 

■J  =  0. 

TO  Tr  4(106 


U  Vi;  K)  U  O  L>  u 


600=5  CON  TIN  LF 

IFfvC.fS.O 1  CO  T»  6307 

*•  F  NO  IT  r  B  ATI  0*i  F  OR  y 


IOC  CON  T I NLr 

U*S  =  I.'  *  CdVE0 

-  ryo  »  rONVFR 

WRITE<f.3101  U.U**S.EyO,EXOMKS 

n  FO«m4TC?X.»XS,JP0  =  *.Fl?.4,3x,.ySI.,F0(',K,El  =  *«Fx?.ft,;y, 
t  »TXSUFO=  *,FI2.*,3X,*TrSU?0«Ml'S)= 

.....  CONFUTE  SIGNIFICANT  VAL"E  V  ..... 


IF  t’J  .S'1.  1.)  It=  E X « 

EVT =  .“*COC 

CALL  NCA-'M  i  (E«T  ,SAM«  IE3  1 

IF<  IER  .ES.12?.0R.I£R  .EO  .HOI  VRI tf(s ,5C 0 II F 5 
'TO  F0R»MT(4fu  rF9C^  IN  GA<*"1  F!JNCTT"N  COMPUTATION.  roo^Ar  ,!») 
AO  s  AvaFA.U 
AUC  =  A  0  *  *  C 

CALL  *  CO  AM  <  ALC.  FVT  .  Pp  .1  F  1 
IF  fC  .LT.  3.)  «  =  l.-CR 
J=3  .»«  l.-PR  l.SAP/CGA»»,A  .*.<*D0A  1 
UVS  =  l  *  CONVEX 
•  pITE<C.3121  U.U-S 

312  F0<»-CT«2*,.ySU?S=  ..Fl2.»,jy,.ySLi:,Rt''l<S>=  «,fi?.4/i 

6*80  CON  T  INL,r 

ySU5F»  -  y  ? l.i 8 0  .  CO‘!VF» 

RSX  «  ?L  )  =  X  SU9G 
95xy(IL»  =  *roep'* 

EXSy  =  Exo,  .  CONNER 
jOITE<6.T3<M  XS0SG.»?CF3«'.FXG*£yvGM 
303  Fov<ATI2X,7HXSU8r-=  ,  F 12  .  4  ,  ?x  ,  1 2H  ysUD  r  f  '<  D  =  .F12.4.SX, 

T  8HT  XS  L8  f=  ,F  1  2  .4 , 3X  «  1  3,;TX  SUE*  G  I v  KS)  =  .F10.41 
C 

c  .....  SOLVE  FOR  ASY"PTOTIC  VALUE  CF  y  USING  VEG  O'r  I N 

C  WFTPOO  S 

r 

FLOG  =  EL/G  A  ,,MA 

IF  (C  .LT.  0.1  rLOG=-EL~r 

US  s  E  >£ A 

s  fOC*  A“9CA)*»C 
NE  =  0 

CALL  NFS  II  ,U. **0) 

C  ..  START  ITERATION  L^O 3  ** 

9010  IF  IU  .LT.  3001  GO  TO  =801 
U  T  =  0  . 

SO  TO  P8  02 

R801  'JT  =  IJ.«'*»Excc-ui 
r  8  02  ,J=C'',’OC4ELCS«UT 

CALL  WEG(2,U»0.NC1 
IF  IU.GE.0.1  GOTO  6011 
'JOaUO*C(  11 
U  a  UO 
J=l  .F-10 
6011  CONTINUE 


NT 


o  o 


•IF  =  NF  *  1 

IF  «*£.LT.6!!1  30  T-  6"1? 

JBITE  (4  ,VS') 

'Ofl'UKI'H  ‘10  CCNVEPOENrE  OBTAINED  F'>o  »CU6FA> 

J  =  0. 

.'■»ITFt£.?0*>  U«  C1C  t  A  WB0  A  .C 
T  3  6H3 
631?  CONTINUE 

TFfNC.EG.ai  00  TO  6M3 

•  •  ENO  TTfoatioN  fob  x  >t 

0  A 

6013  IJ=U«»0CC '»“or4 

u*s  =  i  •  conve  0 

EXS**  s  fXGA  •  COM  WEB 
JPITEt6.631»»  U.UMS.FXG,Fxr a 

631*  fob«at*OX,*X«USOA=  ♦  ,  CJ  J.  A  .  2  X  ,  *X  rUDr-  A  t  “V  c  )  r  *,F1?.*.?X. 

t  ♦TXSUFP»=  ♦*F12.‘»«?F,»TYS'!BOA<'*f<5|=  «,F12.4/) 

6',33  CONTIMLE 

XO'JOOM  =  XOU  =  n  *  CONVE? 

®ox' id  s  ecu3? 

»nx»<Il>  s  YSU5C“ 
rXO“  =  1*0  »  CC*!'/FR 

•IPITEf6.f015)  X'MBO.ySU?r"*.EXD*E  vn« 

5115  format  »2»,*y?UEC=  *  *  F  12  .4  ,2  X  .*  XS'.'B0  (  “*?  )  =  »,F1?.4, 

S  3X,*TXSIJ30=  ..Fl2.4.3X.«TXSUPr,(»k’S)-  *,r12.«) 

C 

c  ««**.  SOLVE  F  CP  ACY*°TOTIC  VALUE  FF  X  'JSINC.  JEOSTBTM 

C  MCTU00  0  *«**« 

C 

•ju  =FL  cwa) 

MO  =  EXOA 
■j  -  ua 
■'E  =  0 

CALL  WEC.a»U»4t"> 

C  ».  STA  °T  ITERATION  L00c  *• 

5043  AU  =  APBCA*U 
11JC  s  AU*«C 

CALL  ’•rr,AY(AUC.'*,PR«  TE> 

TF  «C  .LT.  3.)  PP  =  t.-FO 
U  =  U*LU 
EL3P  =  EL»°R 
UUU  =  f.LPP 

!F< ELPo.LE.3"0. )  UUU = ALOGf E»° ( flpr ) - ; . ) 

?023  U  =  IJ-LUU 

CALL  WEB  ( 2  •U«0»f!C) 

IF  TU.rE.o.j  r-o  TO  604  * 

>ju=u  "  ♦!:  <  i) 

U  =  IJO 

5043  CONTINUE 
NE  =  N  E*  1 

IFANE.LT. 50)  G"  TC  6044 
■JPITE  C5.23PP) 

20 RB  F3RMATA35M  NC  CONVEPPENCE  ORTAINEP  FCP  X'UPCAJ 
M  r  0. 
r.n  TC  60  45 

6044  IFI*IC.f<J.9»  'r  TO  63*0 
C  *•  END  ITE*  ATI  ON  fdd  *  .. 


c  oa 

=  U  ♦  SO'JVEP 
ExOA«  =  EXCA  *  C^NVE® 

S‘4S  *<UTr{  £.«')«<■  1  U.U'SEXCA.EXDA'* 

6346  c0,»"AT*2X,*X':USCA=  *  .  F  1  ?  .4, 2  X  .  *X  SUP  3  !  I  S  )  =  *  .FI  ?  .4.  ->  X  , 

i  »TXSlJ?OA=  *.Fl?.4.?X..TySUB,,'A<'«<S>  =  •»F12.4///J 
10303  roMTIMUE! 

XSITF<6, 60471 

6  34  7  FORMAT  f  1^-1,  r’,4H7If,rf,,T  >  .  T?0  ,=H  TJ«F  fHP  E  >  ,7  79  .  7H‘_  VAL  JE,.*?, 
TSHXSU3C.  167  .FHXSU8G.T77.10HXSUPO  *PKS  >  ,T°  2, 1  0HXSU°G  <<<  S  )  ) 
n!)  9P  1  =  1. M 

esc n>  =  arm  /  &o.o 

xPITE«6.604  8)  ELTt I ) .ESC  1 1 ) .ELLC TO ,PSX  *  T’.°GX<I ) , B ox* f T > ,» Ox “t T 1 

9»  '13*1  T  T*!lr 

6  0  4  8  F09“AT  *1*  ,T;,F13.4,Tl?,Flf>.4,T3  2.F10.4.T4  7,F13.4.T6  2.rl0.*, 
ST7T,F10.4,T9?,F1S.4> 

C 

C  "OTTON  =;  oi_0T  FISTOG°A«*  VITH  PLOT?®  DAC  x AGE 

979  :f*  rs<*i  .lt.d  o,3  n  07 

C  COMPUTE  HIST'i'iRA'*  VALUES 

c 

0IV=1 .3 
ICVT  =  " 

if<  rc<  to  j.ei.i >  niv=ou> 

nP  P?  1=1.11 
tC'IT  r  !  CMT  ♦  1 
HE  Z(IC':T)  =  2E*C(1»/0IV 
33  94  1  s  7.x 

ir< icim.Fo.i)  oiv=ctt> 

33  9J  J  =  3.1C 
!CVT  =  I CNT  ♦  1 
0’  T ( I C *•  T  7  =  0E407  !  )/CI  V 

9»  CO'ITIHUE 

o<*  =  0  .0001 
TPo  =  20 

<13  =  13-K 

?*  Top  s  loo  .  13 

IF( PV( JPO, .LT.P-.ANO.FN ( TPP) .LT.?M.A"3.7* r°Fl .r".0 . 1  1"  TO  8* 

I  F<  IPO  .LT. <101  SO  TO  84 
9*  TP  s  I=P 

TF« T". IT. <191  S3  TO  9? 

'll  *  "l*"  ♦  1 

30  93  I  =  <11.10 

90  ?(I  1  =  0.3 

90  30NTI9LE 

C»LL  I»ITPL9 

C  9°T  10 *1  1?:  NU*  =  E9  OF  PAGES  VEEOED  FOP  pL0T*0FFAULT=1  02UBLF  PGO 
OU9PA3  =  1 

Tc(  ICC  IS  l.NE.J)  *'UKP A5=TC(15> 

**AX  SCA  =  0 
TVAB<1)  s  1«* 

I V*  P  *  2  )  =  *Hr**TST 
IVA»  C3)*EHrGA*J> 

T  VA  9  ( 4  )  =  4«to»t 
I  T*  J 1  )  s  6H3**SE  o  V 
ITP (21  =  6KS  ,5A 


o  o  o  o  n  r> 


I  TP  (  x  ) 

= 

eP'vM»i  AN 

ITP<«) 

6*0  9  A  YL 

7TP  f  5  ) 

= 

£t-f  I*:H  P 

I  TP  (  G) 

z 

©hl^TS 

ITX  <  11 

z 

GH\MLUF 

ITX (21 

z 

till  »AN 

ITX  <51 

z 

f.  w  r  r,  *  v  A 

ITx ( A  1 

z 

6*9 IABLE 

ITX  (*) 

z 

S'** 

ITX  (  G) 

z 

6* 

I  TX  1 1 1 

z 

6P0FNS IT 

rTT(2> 

z 

6H  Y 

ITY(T) 

z 

6P<\Uf*./ 

TTY  C4) 

r 

6 u9  1  N)  /N 

TTY  (51 

r 

SP*?IN-W 

TTY (£) 

= 

6PICTH 

C  tC<5)p?  kSITF.  VALUES  BEFORE  =LOTTTN'*- 

1F<  T-<E1  .F.0.21  J91TE<S*'n9Cltm®(I>,T  =  l.M 
>90  F9B<*AT<lPl,lx,Al»iax»A5,6X.AF,<sX,*4> 

C 

C  OPTtON  Gt  DISTRIBUTIONS  I*  "TS  UNITS 
!*<  ICtEI  .NF.11  9-'  TO  so 
ITX  (51  =  feu  '■»' 

00  SI  I  =  1,19 
X?/*S<I>  s  YN <  I )  *CONV r° 

Tv<  I  j  s  ?« i »  /  cower 
'xisai  s  p’lfi)  /COMvrs 
PN'*S'I>  s  °NfJ>  /CONVfP 

B  1  ’BITE(’)  CXN**Sf  T),2*MI>  ,PN'*SII  ),eN**S  ( »  1  > 

'•-*>  to  i*i 
SO  CONTINUE 

00  91  I  r  1,10 

!F(ICCE>.E0.2>  VBITt<6»3!91)  XN  <  I)  ,  ?( I )  «CN  (  1 1  ,9  f'(  I) 

9191  P0^»7(l"  ,Ar13.<H 
'll  X0ITEf7>  <v*!ft),7<I>,FN<T),9MtT>) 

79  CONTINUE 

CAUL  P  L  9  T  P  9 1 1,9, IV A3 ) 

97  OONT  INLr 

OPTION  T:  Tl»r  EFFECT  PLOTS 

IF<IC{7).N£,1>  On  TC  Q5 

Tr  ST  TC  «FE  7P  T I VE  E^ECT  PLOTS  APE  GPCr 

icGi=a 
I FS  2  so 

TF(NL.EO.I)  GO  T0  95 
DO  99909  1 1 s ? , N L 
IFIRGXItn.NE.PGXtU-ni  IFGisl 
99909  IFIBDX  CT  II  .NF,90»UI-1>  J  TFG2=1 

IcttFGl.n.0.OP.IFG2.EO.P)  SO  TO  °9 
C 

OLL  IMTO(.0 

C  0PTT3*.'  t«:  NUNFE*  nc  PJOfs  NPECED  FOR  PLOT  <  irr  A  (JUT  =  1  DOU°Lr  or-7 
NUXPAC  s  1 


n  n  n  n  n 


IF<  IC<  16).'IE.0>  NUHPAGsI^flG* 

■JPLOT=  2 

>  =  O.C 

C  OP  T 13  N  1  1 :  SCALE  F"R  I  •,OEtrNPENT  VAPIAFLF  C't  FL  OT  ( 'EFA'JLT::  .  1> 
FCA(l)  =  Q.C5 

*F<  IC<  11  ».«P.3»  SCAM  1=10**  (IF  (1  1  1  1 


“AX  SCA 

= 

? 

!  T3  <  1 1 

r 

SHTI-E  E 

I  T3  <  J  1 

r 

6HFFECT 

I  T3  f  .3  > 

r 

6H°L3TS 

TTX  <11 

r 

6HTI'*e'  T 

TTX  <  91 

= 

6~"  POL'3 

TTX  t T> 

r 

6HS 

:ty  ( n 

= 

6-(0Er;I  9N 

ITY  <  ?> 

r 

aH  VALUE 

ITYf  3) 

r 

5W 

ITY  AM 

= 

5H 

I VA  R  f  1  1 

=  9hT  T“E 

:vapi2» 

=  *H»SL,on 

TVAS  1  3  J 

=  £  KXSU9  9 

IVAF  (9) 

r  ah 

C  "’OTTO  !\,  *  ;  Tl'*c  EFFECT  PLOTS  I'l  N*S  UNITS 
C 

c 

IFUCCfl.NE.il  fP  T"  AC 
I VA P  f  2  1  =  Awxr-MKS 
I V  A  F  At)  I  £*YA-J*KS 
ITYCX)  -  tu('t<c, 

SO  61  «,  =  l.NL 

61  JRITE(3)  <ESCtJi«ROXMCJ'*RrxMl J1 > 

0-0  TO  £A 
60  CONTINUE 

00  96  U  =  l.NL 

°6  J 5T  TE I  I)  IESCC Jl.RDX *J>  .DOXC.I)> 

65  CCNTINLE 

CALL  °LCTP»  IT,'  ,IVAR  > 

9=  C  ON  TI NUE 

TF<  TPLCT.EO.C)  G*1  T'  99 

END  loono  L"OF  ♦  •*♦* 


•  *«»«  BLOT  RESULTS 

CALL  3LnTfXLL,0..-3> 

CALL  HISTC 
xN< I P* 1) =F I P  STX 
XN1  IP*?>  =PELX 
F'ICIP*1)=FI9STY 
c«( I=*oi=OFLT 
PNC rP»l)=FI»«TY 
°N( !P*2> iOELT 
CALL  LINECXN.FN.IP.l.O,-') 
JFf  IRAr.E'1.0)  GO  TO  10099 
CALL  LI’lEf  XN.RN.TP.l .5,31 
1C099  CONTINUE 
50  TO  59 

10093  TF( TPLCT.EO.O  STOP 

CALL  PLCTCXLl.e .  .999  > 

STOP 

END 


sueRcumE  oetcpl*  AT.u,r,s*v?AR,3a««  i.^cr •«  r, .  c<> 

CQf*"H  /a  up s»r/v ,i«Br&*c  *r 

C1'*“'M<!?PT/IC°«2C).ICK1 

CC“  vrN  /T,oLlT/i*fl0!;>»ELL<7)»V»£LP.D<l''3)  .  I  ST^e  0 , 1  PL  3T.  I  ®  i  Y  ,  IC  *2  3  ' 
CO'*RC*J  'P°MT.'FORSCIOO  >.FT*<*100>  .OE'lO*  lOO.DE*  T*l<)0> 
t.FT H°<  19C>.3ENT»*13&  >*CC**VE»  .CAL 
rO»»CM/HlST/pIflL'),K.Pff’.«lP;).H«  10")  .Y''£Y,YLEN,YLE»i*C'*VTrS> 

’'IRE*  SIP*;  CO  )  . DEL  "YS«  1  00  > 

DTRFVSIOK  0HEA0*10«> 

°E»L  “ 

OAT  A  0KA0/103*2HO</ 

IF*  IC*23  '  .E  0 .1  )  GO  TO  3 

.....  SOLVE  P"0  "  UEIMO  jPGSTErv  TTPPATTVP  '*ETlnO  ..... 


’IE  =  0 

CALL  U EG  <1 »w ♦ 4.0 ) 

.....  ST£PT  TYPRATICN  LOC3  ..... 


3  r  r3  S  T 1  *“) 

U  a  °S  12  <Y) 

»r*-AT-U/E«*1.5 
coll  arc <2.'*.e,*ic> 

WPIYP<P,AS1)  w 
*M  PCR'UT <1 P0.2»"=»E10.A> 

•IE  a  ME*! 

I FYVC.Lr.SO)  GC  TO  200C 
aoiTP(£,20P) 

206  POR“0T»Y<iH  *!0  C"'*:vERCE'I  CE  ,'ptAI*;pO  F^R  y'UB'’) 
!CK  =  1 
’  E  T 1 J c  N 

2000  IP* «.5E.O.)  GO  f-  2001 
*»  a  .3  00  ! 

20  01  COM  TIN  LE 

I  F*  *IC  ,  c0  .3)  p-0  T"  3 

.....  pno  I  TER  ST  1 3*1  P"°  *  ..... 

.....  COMPUTE  C  ..... 

2  CONTTNL'p 
U 1  =  PEI***) 

'J2  a  PSIK") 

C  a  SORT 'U2/S) 

!c<  T.SE. 0. )  C=-C 
TOC  -  l./C 

.....  CO“PUTE  LSMPP*  ..... 

4  AR30S  a  EXP  f-YBSR.Ul  »OOr) 

5  CONTINUE 
TC<1=0 

IF* IC(23  ).E0.1)  CALL  YLHC^L 
IF* ICK1 .EO. 1 )  RETURN 
URITP»6,101)  CM 


| 


1  n  1  CQR  "AT  <1 1-1  •«  A10  ) 

CXA  r  l.-AL’ 

dRITFI  F,  HGI  CCNVER  .CAL 

ISO  FOR^ATflP  ,3I-CC*:WER=  ,=  !  r  .  *  ,  FX  ,2  C  PC  4  L I  P  F  A  T I  r  »•  e»CT60s  ,=1C.*) 

C  COMPUTc  Xt^RC)  AN 0  DEL»,'vfO> 

00  t?  .1=1. * 

3FL'‘*S(JI)=r)<JI>.C0NVER 
1T  XPKSAJI)  =  X  (JI  )  ♦  C  ONVcs 

C  FIR  D  MK6  V»L"F  Oc  YE  IP 

Y8A  R*K  =  YB  4  P  ♦  K.AL'OtCONVER) 

HRITE<E,*00)  R,*',AL°.  0UA 

300  FOP  ''AT  »//lX,7HN=  , 15 . 7X  .  * PF  =  ,  T5  .  2*  .  7HA  L^A  =  ,'11,*.?'. 

{  2»,0H1-ALPU«=  .013.4) 

VPITc(F.007) 

207  FOP'>aTf/lX.*OEL  ARRAY*) 

JRITEA  c<  29R)  CCPFACI I»«  T.Ott ) »I3  1,»1 
F  OR  ^  A  T  •  1  *•  ,PfA2,I2,lu),r7.4,2X)) 

I  F  (  C  ?‘l  VE  R  «NE  .  1  )  '.RITEfB.OOO 
296  cOR*AT</iX.«rrL  AR  R  A  Y  (“  “’O  )  *  ) 

TF<C'-NVF‘>.NE.1>  WR  ITr  (6 . 29P  )  (  "«= AC f I ) , I ."EL MK5 <  ! > , J  =  1 , R > 
V»ITF(6,301)  YBfiP.YPaR^^.F.T 

*01  =OP,'iT7/iy.fyv9*B=  *  "  12  .4*?Y  *llHYei°  M'VJ  >  =  »?12.P»2X, 

*  *'RS  =  .612.4, 2X.3HT=  .612.4) 

2"3  -0Rv!aTC1Y,2F’4=  ,  61 2.4.2  y,phLAmPCA3  .  r  I  2  .  »  .  *v  ,  me  s  ,61?. 4/) 

z  .....  EVALUATE  njNCTION  ..... 

TF<  IE°  .E-.12S.0  9.IER  .EC  .130)  NRl  Trf  6  .5 0 0  )  t r. c 
F03  c0O“ATt46H  FRR 0°  IN  C4«"A  FUNCTION  COppUYAT I"N.  F30oo=  ,1=) 

I F  (  IC(l=).E0.1.AND.ICI20).NE.l  )  'J°TTF  «  F • 250  >  A«9?A.'*,C 
7  =  3  FOR  M  AT  f  t  *•  ,?  =  HSTACY-''IH7A"  pAR.  A’4?  CAr ,=16 .5 ,44  «r,F!F.5, 

*4H  C=,F16.E) 

COLL  0I0T 
?«0=1.-AL° 

OOC  =  1./C 

CALL  '16A«"A*m,6avMA,  IEc  ) 

AMFS6A'wn?A4CCNYED 
P"*<S  s  R  *  (CONVER..O) 
yft=iof  ftacy/'- 
H00  =°H  I^oA  v 

IFf  IC«20>.E0.1>  «,ET  =  HH>'AXIMU"  LI 
T=< IC<20 ) .F6.1)  HOO=-HK**LI«OCC 
2P7TrfF.?0 2)  M,C.AP59A. A*4**  »R.PMY  =  .yrT »H  nD 
*02  F  0  R  M  A  T  4  /  1X,.*WV-  ,G12  .4.?y,*hc=  ,612,4, 2X,PHLA49D4=  • 

l  ~l  2.4  .2v.l  ?'-LA*PDA  (  M»S  )=  .6l?.4.,*M»s  ,rl?.4,4M0(“F6)r  ,617.4, 
t  /  ,7  H  ,  VIA  ,»n,)«,T»  “EYWCO) 

WORDrS'-PArLEIBH 

«RI  TE(  6,  lit)  >  U3RC,yO®P 

110  FOPNAT  (T.*,PHVALUE  Oc  .T1  9  .RUVALUF  OF  , T*5 . f HMCA  5UR  EO, T* 1 » 

1 A  ft  ,TE6,5HGAMVA,Yft2,8UMrA6UPEC,'reP,AP  .T114.SM6J- 

$./,  T4,  <HRAND0*,,T19,f.HR  A NCO«, T* 5, 9HFR  EQUE  *  CY  ,T51,9HPRE0I  C^EO  •  T66 
i'iHoPE0ICTE0,T82.YHDE*,SI7v,T9R,oMPREr)ICTEn,T114,9HP®ECICYE?,/. 
$TJ,  »PV  AR  IAELF,Ti'’,«HVAR  T8PLE,T51  ,9HF  0  £  CUFI'C  v  .  T66 ,9HF®  EOUENC  Y  , 

8  TP2 , 7H=/PI N /N .T98.7H  PENC ITY ,*114 ,7H0ENS lTY,/«T2C«5Hf'4KS) ) 
‘-•RITE7F,  J20)  <  x  «I)  ,XYKS»I)  ,F69*T  I),FTWo  f  I )  ,^TH  C  l  >  ,"E  v  "«  T  )  , 

4  HEN  TR  < I ) ,0E  NT  < I ),I=1 ,<) 

120  F OR MAT  IT ?,F8. 4, T19,F 4, 4, T37*r12. 4, Y44,r!2. 4, T^*tri2,4, 

6T79,F12.4,T96,F1?.4,T11 1.E12.4) 

0  E  T  U  ®  N 


ooooo  onooooooon 


SUBROUTINE  HISTC 

common  /HisT/eiNL',.Noe i*"-,? to  3 >  ."pd n  "0 >  .that ,xlen *ylen.c''vti c  > 
CC'*.Mr,N/H  1ST2  /F  IRSTX  «  CEL  <  .CI  RSTY.r'ELY 
DIMENSION  »f2C0)»Y(2O0>.RAYIS> 

N2=NCB  INS*  1 
DR0(N2>=Y«AX 

C  OMPLT  £  SOLE  CCR  x-AXIS 
CRD  IS  THE  ARRAY  OF  t4E  1GTHS  Fop  the  FINS 
YLEf.  IS  THE  LENGTH  I  \  INCHES  CF  THE  Y-AXI.S  REFINE"  I ‘I 

a  DATA  STATEMENT  IN  THE  M  A  IN  PRP'-.cif.  TO  =  E  «  INCHES 
N?  IS  THE  NUMBER  OF  POINTS  in  0«n  !»•  c  ♦!  says  Ev E° y  POINT  15 
TO  BE  BLOTTED 

first*  the  first  y  value  annotated  on  thf  Y-axis  is  ret'j»ned 

IN  THT  H2+1  LOCATION  OF  OR E 

DELY  the  SCALE  VALUE  -  AMOUNT  OF  DATA  »fP  INC“-  IS  RETURNED 
° Y  SCALE  IN  the  N2»?  LOCATION  "F  CPC 
CALL  SCALE(ORDtYLEN.N;.*l > 

FIRSTYtORQTN?*! ) 

DEL  Y  =0  PO I N  2*2 ) 

TITLXslOHEXC'JKSIONS 
TITLY  =  thof-;SITY 
ORD(N2>=0. 

SCALES  SET  U=>  BY  AXIS  TO  2E  USED  n  Y  LINE  TO  PLOT 
PnIISTS  -  X-AXIS  LENGTH  IS  CONSTANT  ARC  IS  DEFINED  IN  a 

DATA  ST ATFYEMT  IN  THE  “*IN  PRC?RA«  to  BE  10  I*'CHCS 
THE  AXIS  ANNOTATION  WILL  RE  DONE  WITHOUT  the  USE  OF 
AXIS  SO  THAT  THE  NU“CERS  LINr  UP  WITH  THE  BINS  Dp  AWN 
FOUNT  =  1 
"  =  1 

X(N)  =  3.0 
N  =  N*1 

IF  IBIALC.ro. 0.0  GO  TO  5 
V  = 

X  <N  )  =  PINLD 
■12  =  N2*l 
5  DO  10  I  =  1  ♦  N  2 

XCI  )  =  X(!-l)*RI*.W<  YOUNT) 

FOUNT  =  xOUNT* 1 
10  CONTINL'T 

CALL  SCALE <X «XLEN »N? »*1 ) 

FIRSTX  =  X (  N  2  ♦  l 1 
DELX  =  x  CN2*2 ) 

UNI  T  =  1  . / o E L X 
C  PUT  PEN  AT  ORIGIN 

CALL  PLOT<0.,0.«3> 

JO  20  1=1. NS 
STEP  =  *«I>»LNIT 

C  DRAW  AXIS  PORTION  FRO*  CU°R  EN  T  LOCATION  TO  END  OF  NE*T  BIN 

CALL  PLn  TI  STE°.  0  • .  2  > 

C  ORAN  TICF  BARK 

CALL  PL0T<STEP.-.1»2> 

C  ANNOTATE  AXIS  WITH  NUMERALS 

CALL  NUMBER ISTEP-. 21 .-.21..1.XII >,0.r»?) 

C  RETLRN  PEN  TO  AXID  LINE  WITHOUT  DRAWING  LTNF 

C  CONTINUE  DRAWING  AXIS  AND  A  N*'  OTA  TING  IT  FOR  THE  °E  '*A  I  NnER 

C  OF  THE  VALUES 

20  CONTINUE 


title  the  axis 

C4LL  SYMBOL! A.2S.-.32,.  14, YT  TLX. 0.0.10 
CALL  A YI<<0.,0.,TITLY,*7,rLEN,?3 ..FI RSTY ,DELY) 
RAY<1>=10H*RAYLETGH 
R AY  <  ?)  =  10HCIST. 

RAYt3)=10H-CC'>PUTE0 
»AY(a)  =  10«YALUE 
RAYIS)  r  10H 
X°AGE=XLEM-?.a 
YPA  *E= YLEN-.5 

CALL  S  Y«eCL<''0AGE.Y°AGE..l»,»A  Yt  1  >,0 .0 ,1 0) 
YP43E=XP  4GE*1  .<• 

CALL  SY‘JBOLfXPAPE.YP  AGE  *  .  1 4  ,R  A  Yi  2  >  ,  0  .0.10 

xPAGE=XLEM-2.8 

YPAGE=YLCN-t. 

CALL  SY-eOL<yPAGE.Y°AGE..l».!>AY(  3), 0  .3.10 
XPA3E=»P4G£*1.4 

CALL  SYYP0L<XPA5E»YPAGE..lA.[>AYC«).n  .3.1'’) 

X  PA GE  =  XPAGF*1 .A 

CALL  SYYPOLI  XPAGC, Y=AGE,.1A,RAY(5),C  .0  .1") 
xPA  3E=XLEN-2.H 
YPAG£=YP  AGE-. 5 

CALL  SY“rOLOPAGE.YPAGE..lA.C'*NT,3.P  ,EO 
X  <1  >=B  INLO 

y<d=o. 

X  < 2  )  =6  HLP 
Y«2)=0PC<1> 

J=2 

;n  i  ist.Nipnp 

j=j*i 

XT  =  X  <U-1).PT»IW{I) 

X(J)=XT 

Y<J)=3®0(I) 

J=J*1 

X<J)=XT 

Y<J)=0°C(I»1) 

CONY  INLE 
X  (J*l)  =FIR$TX 
X)J»3)sOFLX 
y<J*1)  sF  1RSTY 
YCJ»2)sOPLY 

CALL  LINECX.Y.J, 1,0,0 
RETURN 


n  o  o  o  o 


PUNCH'")  PSTt'M 

FUNCTION  PSI  EVALUATES  T»E  OIC-AP-A  'UNCTION  RSI»*>  F^R 
AREU'»E*T  *.  PSI«H>  IS  THE  DERIVATIVE  VPT  M  OF  THE  Lt)  OF 
THE  SAo^A  FUNCTION. 


REAL  N 

OAT A  CIS .Cl 126 /.IS 66 666  67, .00701 650R/ 

C 

t  i  =  a. 

'4  -  •" 

IF<  J.GT.tT.)  GO  TO  2 

C  •‘.LF.l? 

"  =  10.  -  V 

N  =  •«  ♦  *' 

0  =  “  -  1 

CO  1  I  -  1,N 

1  Tl  =  T1  ♦  I./(r  ♦  I) 

C  “.ST.1J 

i  5cw  =  i. /vi 

’ey  2  a  pcu*»2 

°SI  a  SLCCM’J  )  ♦  .5*R  CW»  f-1 .  ♦  PC W« ( »C 1 6  «  RCV2*<.1*CIS  -  C  ll?t» 
.  RC'V2J  )l  -  TI 
RETURN 
END 


^UNCTION 


FLV’CTIOK  P$I1  FVAL'IATFS  T mc  TOIr-A>*f*a  FUNCTION  “FT  f«) 
Fee  AF5t'M£,-.lT  M.  PSTIf*>  ts  THF  OE»IVATJVF  U’T  P  ?c  TUF 
OIGAWMA  FUNCTION. 

PFAL  M 

Pi  AT  A  CIS  .C130.C142/.!66S££667,.07337  *3?*  .  .0 2 3S1 9F24/ 

Tt  s  0  . 

VI  r  « 

IPf ..57.13.)  TO  2 
'*.LF  .  13 

'■  =  14  .  -  '4 
J  s  ”  «  " 

0  s  «  •  1 
00  1  Isl.U 

ti  =  ti  *  i./<oi>*«? 

v." T .1 3 

°CW  r  i./vi 
’CJ2  =  SCV*«2 

"Sit  =  ®CV»f l.*9CW*( .«*’CV»tClS  ♦®CV?*(-C1?9  *’CV2*1C1A? 

-°C«2  »C  13  0  )  ) )  > )  *T1 

’FT  UP*) 


FUNCTION  PPI2T*) 


FUNCTION  PSI2  EVALUATES  THE  TE  TP  AGA “H A  FUNCTION  °E T  f"» 

F  CP  »PGUHE*'T  M.  PSI2C“)  IS  THE  ?NO  EEP I  VAT!  VE  VPT  »  "c 
TUE  OT^AHFA  FUUCTICAi. 

PEAL  H 

OAT  A  C16.CE6/.166666EG?,.a3?333333/ 

T 1  =  0. 
u  -  „ 

IF«4.GT.13.1  -0  TO  2 
“.LE.13 

M  =  14. -y 
i  z  <*  ♦  *! 

D  =  M  -  1 
DO  1  1=1  .N 

T1  =  Tl  ♦  2./<0  •*  I  >  «»3 
«.GT.13 

PCV=1./V 
PCV2  =  PCU**2 

PSI2  =  »CV2*»-1.  ♦  RCVI*C  -1.  *PCV*C-.*  ♦PCV2MC1G  -C16 

.  *PCy2*f.3  -C56«  PCW2) ) ) ) ) )  -Tl 

PrT'J°N 

CN3 


no  o  o  o  o  o 


SU8R"UT1>'E  NE*<l.XNPl»J.\> 

UF3STEIN  ITERATIVE  SOLUTION  METHOD  FOP  X 

N*1 


IN  m  AL  I?4  T  ION 

IF# I .NE.l)  TO  2 
J  K  =  1 

x£N  =  XNOI 
XTE MO  i 
X  °  =  X  VP  ! 

RETURN 


2 


FIRST  ITERATION 

IFf * .NE. 1)  SC  TO  4 

*  !F<  APS<Xe-XNPl)-XTEMR*A!sS(X0>>  5,8, P 
a  VP  =  X*'°! 

XbN“l  s  X3N 
XPN  :  X*’P! 

K  s  2 

7  XN  =  X‘,»l 
XN»1  5  Xn*l 
N  s  3 
RETURN 


SUCCEFOI*'3  ITERATIONS 

4  XPNPl  =  f  {X*'°1«»PNIX1  )-(  X-I.XBN)  >/ CX-|P1*XR> '•1-XN-XCN) 

IFf  I  A0S<  XP-XP*IPl>>-<  ^?S  (XTF*,C»XP  >  >>  F.F,6 

a  XP  =  XE'I°1 
x«?N'«l  s  X8N 
X°N  s  XPNPl 
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return 

END 
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RETURN 
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XLOSUMsO . 

X  osuw=o. 

OLNS'JM=0  . 

SU'*0  =p. 
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OTNSRDC  ISSUES  THREE  TYPES  OF  REPORTS 

(1)  OTNSRDC  REPORTS.  A  FORMAL  SERIES  PUBLISHING  INFORMATION  OF 
PERMANENT  TECHNICAL  VALUE,  DESIGNATED  BY  A  SERIAL  REPORT  NUMBER. 

(2)  DEPARTMENTAL  REPORTS,  A  SEMI  FORMAL  SERIES.  RECORDING  INFORMA 
TION  OF  A  PRELIMINARY  OR  TEMPORARY  NATURE.  OR  OF  LIMITED  INTEREST  OR 
SIGNIFICANCE.  CARRYING  A  DEPARTMENTAL  ALPHANUMERIC  IDENTIFICATION 

(3)  TECHNICAL  MEMORANDA,  AN  INFORMAL  SERIES.  USUALLY  INTERNAL 
WORKING  PAPERS  OR  DIRECT  REPORTS  TO  SPONSORS.  NUMBERED  AS  TM  SERIES 
REPORTS;  NOT  FOR  GENERAL  DISTRIBUTION. 


